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1. Introduction. Because noise is an inevitable and undesirable companion of in- 
telligence transmitted or received by electronic systems, it is essential for any proper 
theory of communication to provide suitable methods for studying the physical prop- 
erties of a noise wave and its interaction with a desired signal. On the one hand a suc- 
cessful technique of measurement is required to control or minimize the noise, and on 
the other an adequate theory is necessary to guide experiment and interpret the data. 
Accordingly, the purpose of the present paper is to present a number of new results, 
obtained by the analytical methods developed in recent years, [2-15, 17-20] for the 
following important problems. (In all cases considered here the noise is assumed to 
belong to the fluctuation type characteristic of shot and thermal noise, which are de- 


scribed by a normal random process [5, 10]. Impulsive noise, such as atmospheric and 
solar static, is not treated, although the general methods of analysis remain the same.) 
Our interest here is confined mainly to amplitude-modulated waves, specifically, 

(i) carrier amplitude-modulated by noise: this problem considers the amplitude dis- 
tortion by noise of the carrier wave as the mechanism producing the modulation, and 
the important case of over-modulation is also examined. This example is of particular 
interest when normal random noise is used as an approximate model of speech,f or as 


a form of interference. 

(ii) carrier amplitude-modulated by a signal and noise: this is a frequent case in 
practice where a certain amount of noise accompanies the desired signal in the process 
of modulation. Included also is the problem of speech (normal random noise) accom- 

*Received April 18, 1951. The research reported in this document. was made possible through 
support extended Cruft Laboratory, Harvard University, jointly by the Navy Department (Office of 
Naval Research), the Signal Corps of the U. S. Army, and the U. 8. Air Force, under ONR Contract 
N5ori-76, T. O. I. 

tRecent experiments of W. B. Davenport, Jr. (‘A study of speech probability distributions,” 
Technical Report 148, Research Laboratory of Electronics, (MIT) August 25, 1950) indicate that 
speech is more satisfactorily described statistically in terms of an impulsive or “static” noise model, 
where overlapping among individual (and independent) pulses is assumed to be small, of the order of 
30-50 per cent of the time. This is different from the usual model of fluctuation or norma] random noise, 
which assumes complete and highly multiple overlapping between the elementary transients. However, 
because normal random noise has the great advantage of mathematical simplicity, its use as a speech 
model seems justified on this and physical grounds, at least as a first approximation. 
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panied by noise; the two noise waves are assumed to be uncorrelated. In a later paper* 
is considered 

(iii) s¢multaneous amplitude- and angle-modulation of a carrier by noise: here the 
modulating noise waves are correlated, and there may be in general a phase lag of one 
modulation with respect to the other. The results are of particular interest in connection 
with the problem of the noise in magnetron generators, in which a simultaneous ampli- 
tude- and angle- (i.e., phase- or frequency-) modulation of the oscillations due to the 
inherent or primary noise of the tube is known to occur. 

[A discussion of the problem of carriers angle-modulated by signal and noise has 
been given elsewhere in a recent report (D. Middleton, 1)]. 

We remark further that, apart from the specific applications to problems (i)-(iii), 
the results are needed in the general theory of noise measurement, for here the central 
objective is to be able to determine by measurements on the output wave, following 
various linear and nonlinear operations (such as amplification, rectification, clipping, 
of the original input dis- 


; 


mixing, modulation, discrimination, etc.), the “structure’ 
turbance, i.e., whether or not it is an amplitude- or frequency-modulated wave, how the 
noise and signal occur together, and other qualitative and quantitative data. 

The quantities of chief physical interest are (a), the mean or steady component of 
the disturbance, (b), the mean intensity of the wave, and (c), the spectral distribution 
W(f) of the mean intensity. This latter quantity is in fact sufficient to give us the other 
two; the mean intensity is obtained as the area under the spectral density curve W(f), 
while. the (square of the) steady component is given by the constant (or frequency- 
independent) term in the expression for the spectrum. It is assumed that we are dealing 
with a stationary (and ergodic) random process, namely, a process for which the under- 
lying mechanism does not change with time. Then time averages and ensemble or 
statistical averages are equivalent, [20] to within a set of random functions of prob- 
ability zero, so that if we represent our stochastic, time-dependent disturbance by y(é) 
we may write the steady component: (a) as** 


l 


T ) 
(y(to))av. = lim T™ [ y(to) dto = (y(to))..av. = [ yW,(y) dy, (1.1) 
T-@ “0 J—@ 


and the mean intensity (6) 
‘ 
(y(to)*)av. = lim T™' [ y(to)” dtp = (y(to)”)s..v. = y W,(y) dy; (1.2) 
Ta “0 « 


W,(y) dy is the probability that (at any initial time ¢)) y lies in the range y, y + dy" 


The moment of greatest interest, however is given by 


7 
(y(to)y(to + 2))w. = lim ia [ y(to)y(t, + t) dt . (1.3) 


T- 2 “0 


R(t) 


The quantity R(t) is the auto-correlation function of y and may be found statistically 
y, ; t) is known; W, has the following 


when the second-order probability density W.(y, , 

interpretation: 
*We distinguish here between time and statistical averages by { )sv. and ( ), av. respectively. 
** Accepted for publication in the Quarterly of Applied Mathematics. 

The equivalence of these averages follows from the ergodic theorem [2, 20]. 
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W(y. , Yo 3 t) dy, dy. = the joint probability that at some (initial) time t) , y(= y,) lies 
in the range (y, , y, + dy,) and at a later time ¢, + t, y(=yz2) 
falls in the interval (y2 , y2 + dye). (1.4) 


Because time and ensemble averages are equivalent here, Eq. (1.3) becomes 


R(t) = (YrYo)e.ar. = [| YY2Wly. , Yr 3) dy: dy, (1.5) 


JJ —@ 


and since the process is stationary, the initial times ¢, do not enter: one is concerned only 
with the time intervals (t) between observations. 

Knowledge of the correlation function R(t) is important, for by the theorem of 
Wiener [14] and Khintchine [15] the mean intensity spectrum follows at once as the 
cosine Fourier transform of R(t), namely 


Wi =4] RW cosatdt, (w= If) (1.6a) 


with the inverse relation 
R(t) = | W(J) cos wt df. (1.6b) 


To determine the desired energy spectrum W(f) the simplest procedure is first to obtain 
the correlation function and then apply (1.6a). Note from (1.6b) that setting ¢ = 0 in 
R(t) gives the mean total intensity of the random wave, namely 


ro) 


RO) = | W(f)df = lim [| Yiy2Wolyr , Yo 5 t) dys dy2 = (Yy’)s.v. (1.7) 
0 t-0 J /-—@ 


which is the area under the spectral distribution curve W(f), as expected. On the other 
hand, allowing ¢ to become infinite in R(t) yields the steady component (y),.... , since 
lim,.. Ws(y, , yo 3 0 = Wil(y,)Wi(y2), so that (1.5) becomes 


lim R(t) = ff nueWin)Wilys) dy dys = Ww. » (1.8) 


from (1.1). For a pure noise wave (y),.,., Vanishes, as there are no steady components. 
However, when y does not represent a purely stochastic variable, but contains steady 
and periodic terms as well, lim,.... R(é) will not die down in time, but will oscillate in- 
definitely. If R(t) is then expanded in a Fourier series, the coefficient of each periodic 
component represents the mean power (or intensity, as the case may be)* associated 
with that component; setting ¢ = 0 in R(é) still gives the mean total power in the wave. 

In a similar way we may find the correlation function for a general function g(y) 
of the random variable y. By definition (3, 10) we have 


R(t) = (gly(to)lgly(to + 2)])av. 


(g(r) 9(Y2))s-av. 
(1.9) 


aw 


[P ouroud Wats 5 ve 59 dyn dy « 





*The mean intensity may be expressed in units of power or mean square amplitude, appropriate 


to the problem in question. 
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The spectrum follows from (1.6a). For problems (i)-(iii) the modulated wave V(t) is 
expressed as a function of a statistical variable y, and the choice of g(y) is based on the 
pertinent physical model which describes the problem. In general, g(y) is not a linear 
function of y, and so the evaluation of the auto-correlation function becomes difficult. 
These remarks are illustrated in the following sections. 

The main results of the analysis of a carrier wave amplitude-modulated by normal 
random noise, or a signal and noise, show that the amplitude-distortion characteristic 
of over-modulation spreads the spectrum but not significantly; the additional noise 
components are due to (n X n) noise modulation products. Furthermore, when a modu- 
lating signal accompanies the noise, distortion of the signal also occurs, and (s X n) 
as well as (n X n) noise harmonics are produced. Expressions for the mean total power, 
the mean continuum power, the mean carrier power, and the mean power in the discrete 
portions of the spectrum are given, along with a detailed treatment of the spectral 
distribution of the wave’s energy. In the following sections, a discussion of the limiting 
cases of weak noise, strong signal, etc., is included, and a number of figures illustrate 
the principal results. 

2. Carrier amplitude-modulated by noise. We represent the /F (or RF) wave by a 
complex disturbance 


g(y) = V(t) = A(t) exp (twot), (wo = 2rfo), (2.1) 
where A,(¢) is a real quantity. The amplitude modulation is specifically 
A,(t) = Ao(l + kVy(d), y=kVy() => -1 


kVx() < -1 


= 0, y 


in which V,y(é) is a normal random noise voltage (or current) and k is a modulation 
index, with dimensions (volts)~*. When the instantaneous amplitude V y(t) is less than 
—1/k, over-modulation occurs, and the signal generator does not oscillate until V y(t) 
is once more greater than —1/k. Since we assume a purely normal random noise, large 
and even infinite amplitudes are possible, and consequently we may expect over-modu- 
lation for a noticeable part of the time, unless the modulating noise is weak. The analysis 
of this and succeeding sections assumes the common type of modulation in which the 
instantaneous amplitude of an oscillator’s output is modified according to some signal 
or other low-frequency disturbance applied to a suitable control grid. Frequently, how- 
ever, a modulated output is produced by applying the sum of the separately generated 
oscillations and the modulation to the input of a (half-wave linear) rectifier. The tube 
acts now as a mixing device, which yields a suitably modulated output carrier wave 
only if the original carrier oscillations are very intense relative to the modulation. Otherwise 
one obtains serious distortion due to the significant additional harmonies generated in 
the nonlinear mixing of the signal (noise), carrier, and background noise. Thus, if a 
mixer is used, Eqs. (2.1) and (2.2) apply here approximately, provided the modulation 
is weak, while (2.1) and (2.2) are valid models for all degrees of carrier and modulation 
strengths when the alternative system of a modulated oscillator is employed. 

Let us consider first the simpler and less general situation in which the r-m-s noise 
amplitude (Vy (é)*).”2, is very much less than 1/k 2'”’; this means that for an overwhelm- 
ingly large percentage of the time the instantaneous amplitude kV y(t) is less than unity 
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and therefore over-modulation is for all practical purposes ignorable. The exact ex- 
pression (2.1) then becomes simply 

A,(t) = Ao(1 + kVy(t)) = Ado(1 + y), (-—2 <y=kVy(t) <@), (2.3) 
provided (0 < k(2p)'”* « 1), where by y we abbreviate (V y(¢)’),.... . To see how strict 
the condition (0 < k(2y)’* <« 1) must be, we ask what fraction of the time a 
(0 < a < 0.5), V(é) exceeds an amplitude V, = 1/k. Since V(t) is normally distributed, 
we have at once 


a 2, » 1 1/2 
= ey Iocan CRP CK V/2W AV = 51 — OC1/KQW)")), 








(2.4) 
(y = (Vy (t)*)..0%.); 
where @ is the familiar error function 
es eae my.) — 4” exp (—2*/2) 
O(x) = 2x | exp (—z) dz, and ¢” (x) = dx” Qn? (2.4a) 


With the help of tables of @ we readily determine k(2y)'”’ corresponding to a chosen 
value of a. We see for example that when k(2y)'”? < 0.6, over-modulation occurs less 
than 1 per cent of the time, and so for most purposes Eq. (2.3) may replace the more 
general relation (2.2) when (Vy(t))..av. < 0.18/k. 

The autocorrelation function of the modulated wave is therefore from (1.3) and (2.3) 


R(t) = (A5/2) Ref exp (—two [1 + k*( Vn (to) V(to + t))av.]} 
(2.5) 


= (Aj/2) cos wt{1 + k*Ro(é)y}, 


in which R,({)y = (Vy(to)Vwllo + 0)s.av. iS the auto-correlation function of the modu- 
lating noise wave. The mean intensity spectrum is from (1.6) 


W(f) = (43/2) 8(f — fo) + AdK? | Roly cos (wo — w)t dt, (2.6) 


where we neglect the contribution from cos (@) + w)t in (2.6), since the spectral width 
(—w,) of the noise is much less than the carrier frequency (= f,). The carrier power is 
unchanged, viz., W,, = A>/2, while the amount of power W, in the continuum is dis- 
tributed symmetrically about f, , with a total intensity Ajk*y/2. [Unlike frequency- or 
phase-modulation (ef. sees. 2, 3 of ref. [1]) there is a larger amount of energy in the 
wave after modulation than before, and all of the additional power appears in the con- 
tinuous spectrum.] Now Ro(é)y cos wet in (2.6) represents the correlation function of 
the continuous part of the output spectrum, which is the same as the correlation function 
that one obtains for a narrow-band of noise centered about f, , and consequently yields 
the same intensity spectrum. There is, however, an important difference between the two. 
In the former the component at a frequency fy + f’ is correlated with the component at 
the image frequency f, — f’, while in the latter there is no such coherence between pairs 
of harmonics symmetrically located about fy . Nevertheless, identical forms of correlation 
function and spectrum occur in either instance because, by definition, these quantities are 
squares of a modulus, from which all phase factors are excluded. On the other hand, a 
Fourier analysis of the noise-modulated carrier (2.1), (2.3) and the equivalent (in power) 
narrow-band noise centered about f) show at once coherence in the former and none in 
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the latter. This can be observed directly on a cathode-ray oscilloscope: for a noise-modu- 
lated carrier with ignorable over-modulation, the instantaneous envelope will vary 
randomly but there will be no change in the phase of the carrier, while in the case of the 
noise band the phase of the carrier will change (relatively slowly) in a random way. 

We return now to the general case (2.2), which includes over-modulation. To represent 
the discontinuities in A,(¢), ef. (2.2), we express A,(¢) in terms of its Fourier transform 


2° dz exp (iz[1 kVy(d))), (2.7) 


JC 


L(t) = —Ag(2m)"' | 
where C is a contour extending along the real axis from — © to +o and is indented 
downward in an infinitesimal semicircle about the singularity at z = 0. The correlation 


function is now 


R(t) = (Aj/2) Red exp (—twot)(4m°)' | 2 * dz exp (iz[1 + kV y(t) ]) 


* dé exp (—7€[1 + kVy(to + op} 


ji 


since no coherence between carrier and modulation is assumed. Here C* is the contour 


wtr 


. stat.av. 


conjugate to C, extending from + © to — © and is indented upward in an infinitesimal 
semicircle about the point § = 0. Furthermore,.inasmuch as V y(t) is real, A,(¢) is also 


= —§ 


A,(t) and we can replace the integral over C* by one in C, setting 7 
The ensemble average in (2.8) may be effected in a straightforward way if we note that 
Vy(lo) (=V,) and Vy(to + 8 (=V-.) are normal random variables, their joint 


vil 


and so A,(t)* 


since 
distribution is given by a relation of the form 


WAV, , V2; 0) = [2e¥(1 — 15)']"' exp {-—[Vi + Vi — 2V,V.r0]/2(1 — r3)¥}, (2.9) 


end 7; = (V.V2),.<0 JV ees y(t)/W is the normalized auto-correlation function of 
the modulating noise, whose mean intensity spectrum is W(f)y . Substituting (2.9) into 
(2.8) and observing that the statistical average yields the characteristic function (ef. 
Eq. (2.16) of ref. [10]) for the noise F(z, £; Jy = exp |— 3hk°y(0)(z* + £) — zéy(d)}, we 


have finally 


R(t) = (A5/2) Re§ exp (—twot)(4x°)™ I. z ° dz exp (iz — k' yz" /2) 


| & ? dé exp (i — k’y#’/2 — woah 


JC 


n 


(A2h? ,/2) cos wot + (A2/2) >! : Vrol i h2_, COS wot, (2.10) 
n=1 ° 


where 


ho.» = (2n)~? | 2"? exp (iz — k’y2"/2) de = —i "(Ry)" 


fofa-it 1 \//-x 
4aFif 2 °2 sig) /™\ 2 ) 6.1%) 


2 + | 5 ee )/ (2=)) 
+ (2k* yp)" 2 if (2; 9’ Qk? yp I 2 f 


/26>(n—-3)/2 
oie 
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from Eq. (A3.17) of reference (10); ,F, is a confluent hypergeometric function. Specific- 
ally (ef. (A.9) of ref. 10), we have for the amplitude functions ho,,, , 


hoo = [1 + O([2k?y]~'”*)]/2 + (k?y/2m)'” exp (—1/2k’ yp) (2.12a) 


ho, = a[1 + O([2k’y]~'””)]/2, (2.12b) 


how = (IYI 9" (VT), (w= 2,3, 4,5, ++) (2.120) 


(the definitions of 0 and ¢@™ are given in Eq. (2.4a); see also Appendix III of ref. [10]). 
By Eq. (1.6a) the mean intensity spectrum is the Fourier transform of (2.10), namely, 


W(f) = (Atho.o/2) 5(f — fo) + Ack?(—hi..)¥ [ To(t) cos (w. — w)t dt 


(2.13) 
+ Aok*y > ue. ; To(t)” cos (wo — w)t dt. 
n=2 : Jo 
When the modulating noise has a gaussian spectrum, Wy(f) = W, exp (—w’/ws), this 
becomes explicitly 
W(f) = (Abho,0/2) 6(f — fo) + 2? Ack? Por {ie exp [—(wo — w)*/ws] 
(2.14) 





o n-2 2 —1/2\2 
eo (Tk? yy”) 2 
+ iin? exP [=o — «2)*/ne] f. 
Figure (2.1) shows typical intensity spectra for a number of values of (2k*p)'” between 
0 and o. We distinguish two limiting cases: (2k*p)'”? > and (2k*p)'”? — 0; the auto- 
correlation function (2.10) is accordingly 


9 2 > Ask’ y us 
(2k*y ~~): R(t) = 008 wot-41 + 9 Tol?) 





(2.15a) 
“ —serg(t)”™**(2n)! is 
+ d 2’"n!’(2n + 2)(2n + 5} + O([2k"y)"”), 
” AEN cos wut-{ra(0] 5 + sin™' re 
(2.15b) 


+(1- ‘} + 0([2k?y]'””) 


for the former, while in the latter instance we obtain as expected Eq. (2.5), with a cor- 
rection term o(e **’*nr,(t)*). 

Equations (2.15) apply when the maximum amount of overmodulation, namely 50 
per cent of the time, occurs, whereas (2.5) yields the correlation function when essentially 
no overmodulation takes place. Additional correction terms may be found in straight- 
forward but tedious fashion. We note from Fig. (2.1) that the mean intensity spectrum 
is here more widely distributed about f, than for the case of ignorable over-modulation. 
The additional harmonics are (carrier X noise) noise products, stemming from the 
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Fic, 2.1. Mean intensity spectrum of a carrier amplitude-modulated by noise. 
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clipping process inherent in overmodulation. However, [unlike the examples of fre- 
quency- or phase-modulation discussed earlier in ref. 1] there is a limit to the spread of 
the spectrum, determined by the fact that over-modulation can occur at the maximum 
but 50 per cent of the time. On the other hand, since the amount of clipping inherent in 
the weak modulation cases is ignorable, no significant spectral spread is obtained, (see 
Fig. 2.1). In any case the spectral spread is relatively small. 










$f = 
We * MEAN TOTAL POWER 
Wig * MEAN CARRIER POWER 
— W. = MEAN CONTINUUM POWER 
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Tt 'y 


je 
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Fic. 2.2. Mean power in a carrier amplitude-modulated by noise. 


Whereas the spectrum requires a series development, cf. (2.13), the total mean power 
W, and the total intensity W. of the continuous part of the disturbance are easily ob- 
tained in precise, closed form from Eq. (2.2) if we remember that kV y(t) = y is normally 
distributed with the first-order probability density W,(y) = (2rk*y)~'” exp [—y?/2k’y]. 
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The power in the carrier after modulation is (ef. (1.1) and (1.2)) 


2 


W,, = (A3/2)| | (+ y)Wily) dy | = (A2/2)h3 0 
se (2.16) 


= (Ao/2){[1 + @([2k*y}'*))/2 — kyo (Tk y]”)"} 


W, = RO) = (A3/2) | (1 + y)?Wily) dy 
v-1 





(2.17) 
= (43/2| 24) ouaevr) + 0] - ever’) 
and so 
W. = W, — Wy, = (A3/2)[1/4 + (k’y/2)[0{(2k?y) 7} + 1] 
(2.18) 


— {O[(2k’y)~"”?] — k’yo([k?y]*””)}?]. 


Figure (2.2) illustrates the behavior of the mean intensities W, , W,, , W. for different 
degrees of over-modulation. (See also Figs. 3.1 and 3.2 when w = 0.) Their limiting 
conditions are instructive: As the intensity of the modulating noise is increased, a 
correspondingly greater proportion of the modulated wave’s power is distributed in the 
(noise X noise) noise sidebands, generated in the process of modulation and the clipping 
due to possible over-modulation, as shown in Fig. 2.2. We remark that the amount of 
power in the carrier component and in the continuum is quite independent of the partic- 
ular spectral distribution of the original (normal random) noise, and depends only on 
the clipping level at which over-modulation occurs, since power in a given band of 
frequencies is proportional to the integrated intensity of the band. One can consider 
also more complicated modulations, such as square-law modulation, viz: Aj(f) = 
{1 + kV,Jt)]’: the treatment is identical with that of the linear case examined here, 
following an appropriate modification of the transform relation (2.7). 

3. Carrier amplitude-modulated by a signal and noise. The results of the previous 
section may be generalized to include the important case of modulation of a carrier, 
Ao exp (tw), by a signal with an accompanying noise disturbance. The modulated 


carrier (2.1) may now be written 


V(t) = Ao(t) exp (twol) 
= Agfl + kV y(t) + uVs(6] exp Got),  kVy+uV, > -1 (3.1) 
= QO, kVy t+ UV s 5 =] 


which includes possible over-modulation when the signal and noise are properly phased 
and sufficiently intense. As before, A,(t) is represented by a suitable contour integral 
[ef. (2.7)]. Remembering that A,(é) is a real quantity, we apply Eq. (3.1) and the results 
of section 2 to the general relation for the auto-correlation function of the modulated 


wave and obtain finally 
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R(t) = (] 2) Ref V(t) V(to + ae 


= (Aj/2) Relexp (—twot)(4x°)™ [2 exp (iz) dz (3.2) 


, [. £* exp (i) F.(2, &; OP 2, &; ts ax, 
where F(z, £; {)y is the characteristic function for the accompanying noise, and 


F(z, g; t)s = [exp (iu V s(to)z + ip V s(to + OE) le cat.eve 
To (3.3) 


=T7;' | exp GuVslt)e + tuVs(to + 08) dt 


is the characteristic function for the signal, with 7’, the period of the modulation. Ex- 
panding the exponents in (3.3) in a double Fourier series and averaging over the period 


T,(=2r w,) we can write the signal’s characteristic function in the general form 
F(z, 3s = > (—1)"enBn(2)B,(€) cos mat. (3.4) 
m=0 


For signals which are entirely stochastic, we may replace the time-average in (3.3) by 
its equivalent ensemble average, since we are assuming throughout stationary (ergodic) 
processes. An example of the latter type is provided when the modulating signal is a 
(normal) random noise, uncorrelated with the background interference, in which case 
F(z, £; ts = exp {—n’¥(0) (2 + &) — w’W(d) s2k}, 
(3.5) 
W()s = W0)sro(d)s , (¥(0)s = Ws). 


Here y(t), is the auto-correlation function of V(t)s . The correlation function of the 
modulated wave is then given by (2.10), provided we replace k*yro(t) therein by 
kvyro()y + wbsro(t)s and k*y in ho, , Eq. (2.11), by k’by + ups . The resulting in- 
tensity spectrum is a superposition of (n X n) noise components, the carrier being the 
only periodic term. 

Let us consider now the more complex situation involving a periodic signal. The 
auto-correlation function (3.2) becomes with the help of section 2 and (3.4) 


R(t) = (A2/2) p » €n(—1)"[cos (w + mw,)t + cos @o — mw,)t] 





(3.6) 
1f 2 2 (—1)"(k’yro(t))” 52 } 
“ oe + x n! ect 
where the amplitude functions h,,,,, are 
ha. = (2x) [. 2"~°B,,(z) exp [iz — k’yz2"/2] dz. (3.7) 


The mean intensity spectrum follows at once in the usual way with the aid of the theorem 
of Wiener and Khintchine, ef. (1.6). Because of the distortion inherent in over-modula- 
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tion we expect that the original signal will be deformed, and because of the accompanying 
background noise we may further predict that not only will there be cross-modulation 
between the components of the noise, but between them and the signal harmonics, 
modified in the process of over-modulation. Accordingly, we observe from (3.6) that 
the term in the correlation function for which (m = 1, n = 0) corresponds to the carrier 


component, and the harmonics for which (m > 1, n = 0) represent the (s X s) signal 
cross-terms generated in the course of over-modulation. On the other hand, the com- 
ponents for which (m = 0, n > 1) are attributable to (n X n) noise harmonics, while 


for (m > 1, n > 1) one has as the noise contribution (s & n) noise products. The mean 
power associated with the carrier, signal, and continuum are obtained from (3.6) on 
setting { = © for the periodic components and ¢ = O for the stochastic part of the 
modulated wave, according to section 1. Note that again, cf. section 2, the power content 
of the disturbance does not depend on the spectral distribution of the noise and signal 
modulations. (See also Appendix ITI of ref. [10].) 

We can determine W, by a more direct method than expansion in a double series, 
which involves one less infinite development, unlike (3.6). The procedure is based on 
the observation that y = kV y(t) and z = uV s(t) can be treated as independent random 
variables, whose first-order probability densities W,(y) and w,(z) are easily determined. 
Thus, for the background noise, W,(y) is a gauss distribution density for which (7°)... = 
k’p, (y)s.av. = 0, while for a periodic signal 


lx 


* »2 


4,(z) = | (2r)~* exp (ize) dé | (2r)~* exp [iuEV 5 (0, )] do. (3.8) 
The phase ¢ is a purely random quantity, distributed uniformly between 0 and 2z 
with a probability density 1/27. Therefore V; is also a random variable, corresponding 
physically to the fact that we agree now to observe the periodic wave at (independent) 
random times. (In this fashion any periodic disturbance can be “‘randomized”’ with 
respect to the observer.) The mean-square value of the modulated carrier (3.1) becomes 


W. = (| V(to) |s.av./2 = (A0/2)(1 + y + 2)”),.av. , (for allw = y+2 > —1), 
oc (3.9) 
= (42/2) | W{?(w)(1 + w)* dw = RO), 
J-1 
where the frequency-function W,"’(w) for w = y + zis given by 
Wi? (w) = | Wi (y)w,(w — y) dy 
(3.10) 


.29 


= | (2r)° exp (—iwt — k pe /2) dé | (27) exp (inEV s(0, )) do, 


since the Jacobian | d(y, z)/d(y, w) | of the transformation is unity. In a similar way 


one obtains the mean power in the carrier, which is 


W,, = (AG/2)ho.0 = (| V(t) |)2.av./2 = (Ao/2) | (1 + y + Zev. |’, 
or é = ee = 
for all u y¥+22-~! (3.11) 
= (Aé/2)4 | Wi" (w)(1 + w) dw , 
J-1 
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since w is real. Note that this procedure does not give the mean power W,,,. in the 
periodic components of the wave, but only the steady or average value of the envelope. 
To calculate W,.,. one must sum the series 


Oe te 1)"hin.0 . 


m=0 
In the specific case of a sinusoidal signal, 
Vs = V> cos (wt + ¥), 


we find from (3.3) and (3.4) that now B,,(z) = J,,(uVoz), and so the amplitude functions 


> 


3.7) are explicitly* 


h, = (2r)™' I. 2” "J Au V2) exp (zz — k? pz" /2) dz 
4 o” k* a , r , 1/2 ym 
2m! ( 3H) ‘[wVo/(2k y) | (3.12a) 


(2), Fi(lg + m + n — 1]/2;m + 1; — u Vo/2k*y) 


qk’ y)’?T((3 — m — n — q]/2) 


Ol 





—<""" (k*y\"""” ee 2472 o2.2.2\€ 
= o (EY) [a Vo/(2k*p)'?]” De emma? Vo/2k*y)*, (3.12b) 
where 
(F ((2¢g + m+n — 1]/2; 1/2; — 1 ‘2k? p) 
a i ale ! ao bon A SO MS) Bs Bia ae. f 
—_ arg + m) 1) r((3 — 2g — m — nj/2) 


(3.12c) 





aye Fi([2g + m + n)/2; 3/2; — 1 ao) 
+ (2/k' yp) r'((2 — 2qg — m — nj/2) a 


The first result, (3.12a), is convenient when the signal is large relative to the noise 
[uV; > 2k*~], for then we may use the asymptotic form of the confluent hypergeometric 


function 


iF,(a; 8; — z) = re) r at + a(a — 8 + 1) 


r(B — a) x1! 


(3.13) 


4 ea + Da- B+ I@-B+2, sof Re(z) > 0. 
x2! 

On the other hand, the series (3.12b) is particularly useful for weak signals [uV¢ « 2k*y]. 

From either expression we can easily obtain the important limiting cases of over- 

modulation and weak noise (2k*y — 0), the latter by (3.13). 


*The details of the integration are given in section 2 and in Appendix III of ref. [10]. See also Ap- 
pendix IV of ref. [3]. 
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Fig. 3.1. Mean carrier power of a carrier amplitude-modulated by noise and a sinusotdal signal 


o. 


The mean power in the carrier is found directly from (3.12) to be in the sinusoidal 


ease Ath,.,/2, which is illustrated in Fig. (3.1) for a variety of values of (2k°y)'~ and 
u’V> . The mean total power W, is in the present instance more easily found from (3.8)- 
(3.10) than from the multiple series expression, where now V'.(0, ¢) = Vo» cos (y + ¢). 
We have 


W, = (Aj/2) | (1 + w)* dw | a Jo(uVot) cos wt exp (—k yét’/2) dt = RO) 


(3.14a) 


= (43/2))(5 + SP 4 HP) + octan yy) — hve 10 Y) 


‘3.14b) 


+ yD (Se) eo OY. 





1952] DISTRIBUTION OF ENERGY IN NOISE-AND SIGNAL-MODULATED WAVEST _ 351 
W, is shown in Fig. (3.2) for representative values of u’V> and (2k*y)'””. The power 
associated with the carrier part of the modulated wave is obtained alternatively from 


(3.11) in a similar manner; we have finally 


W,, = (Ao 24 (1 + Of(2k?y)'7}]/2 — bye [Py] 
(3.15) 


2 »\1/2 - *Vo = (2n-2) 2 ,\-1/2 . 
+ (ky > (473) 7p? (hey) 


which is equivalent to our result A¢h5 9/2. Typical curves for W,, as a function of (2k*y)'” 


are shown in Fig. (3.1). Since W,,.. = W,, + Wz.) one easily finds the mean power 
W 2s) associated with the signal components once W,,., and W,, have been calculated. 
Furthermore, because W, = W, — W,,., , we obtain also the mean power in the con- 
tinuum, without having to sum the doubly-infinite series of the direct expansion (3.6). 

lor a given amount of signal (uV, fixed) the power in the periodic and in the carrier 
components of the modulated wave increases with the amount of modulating noise, as 
shown (for W,, only) in Fig. (3.1). The energy available in the signal becomes inde- 
pendent of the noise; however, the noise is then relatively so great that the signal is 
quite ignorable. This is easily seen from (3.12) in the case of the sinusoidal signal (—hj,o) 
when 2k°y —o. Depending on the strength of the signal and noise relative to the 
amplitude A, of the unmodulated carrier, some of the remaining signal power is dis- 
tributed in (s X s) “diserete,”’ periodic terms (m > 2, n = 0), which represent a dis- 
tortion of the original sinusoid, attributable to the clipping inherent in over-modulation. 
l‘urthermore, as the noise becomes more intense, over-modulation occurs a significant 
fraction of the time, up to a maximum of 50 per cent. The amount of signal power in 
the modulated wave is then one-half that in the original modulating signal, since on 
the average the comparatively weak signal “rides”? on the stronger noise half the time. 
Additional noise is also generated by the clipping of the wave, and these new noise 
components appear as (s X n), (m,n > 1), and (n X n), (m = 0, n > 1) terms, pro- 
duced by the cross-modulation of the signal and noise harmonies. Only the (n X n) 
terms are important when the noise is strong compared to the signal. In a similar way 
one finds that for weak noise and strong signals, the (s X n) noise products are sig- 
nificant provided ywV, is noticeably greater than unity (over-modulation of the carrier, 
due essentially to the signal). If uV, is less than unity and there is little noise, we expect 
a negligible amount of over-modulation, and the original signal consequently suffers no 
appreciable distortion. The same argument applies to describe the variation of the 
mean-total power W, and the mean continuum power W, = W, — W,,., with different 
amounts of signal and noise modulation. Again, when the noise is the dominant factor, 
both W, and W. become indefinitely large as 2k°y ©, ef. Fig. (3.2). We note also 
that W),., approaches a fixed limit, independent of the amount of modulating noise, 
which is one-half the original power in the carrier and modulating signal. Most of the 
wave’s energy goes now into the noise continuum (s X n,n X n), as a result of the 
very heavy (—50 per cent) over-modulation. When the noise becomes weaker, corre- 
spondingly less of the modulated carrier’s power is distributed in the continuum. In 
general, different degrees of modulation (i.e., different values of « and k) will, as ex- 
pected, critically affect the magnitudes of total, periodic, and continuum powers. 

The explicit calculation of spectra, based as it must be on Eq. (3.6), is far more 
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tedious than the determination of power. A sketch illustrating typical spectra is given 
in Fig. (3.3); the precise calculations are reserved for a later paper. 

With the help of the properties of the hypergeometric function one can obtain other 
interesting limiting expressions for the power and the spectrum.* For example, when 
there is ignorable over-modulation, (u7(V §)..av. + k°(V n)s.av. K 1), Ann» (mM > 2,n > 1), 
approaches zero, and one has as expected 
R(t) = (A2/2) 11 + rl + eV s(to)Vis(to + D)av.} COS wo, 

(3.16) 
(u(Vs)e.av. + kV n)s.av. K 1). 


For no signal at all the results of the preceding section can be applied, while at the other 
extreme of no modulating noise (¥ = 0) one finds R(é) from (3.6) and (3.7) on setting 
0’ In the specific case of sinusoidal modulation the amplitude functions are** now 


yn = 
«6 = | J (Voz) exp (tz) dz/2nz 
Jc 
_ eer, : oF ({m — 1]/2, [—m — 1]/2; 1/2; 1/u° Vo) 
7 | r((3 + m]/2)r({3 — m]/2) 
(3.17a) 
2 .F,(m/2, — m/2;3/2;1 ae} , 
“ie : 77) ee a 1 < pV, 
+ uVy, {2 + m]/2)T((2 — m]/2) (1S #Vo) 
= ==—-] , m = @, 
ip V,/2, m = 1, 7. (0 < pV, < 1). (3.17b) 
Q, m > | 


The mean total power is found from (3.14a) to be, (k°y = 0), 


fi vet 2,8 _— wV, >I 
i-——CCO™ 
I 


> l o =] . . 2 Vo ‘ efo q § 
= (A3/2) (3 + * sin z0)(1 + He ) +7 (1 — Vote) - zai} (3.18) 


while the mean power in the carrier is 


W. = R(O) = (A;/2z) 
0O<yuV, <1 


, 


- 
1 


W,, = (42/2) | | (uVoz + 1) 


; pC 
ee 2)? 


> 


2, ™( , 


l as V en , . 
= (A? ae +-sin'z +°°(1-2 al . (2k?y — 0). (3.19) 
The difference W, — W,, now represents the mean power in the continuum, in this case 


*See section 4 of ref. [10]; in particular the case of the biased, v th—law rectifier. 
**See Eqs. (4.19), (4.20), and (4.22) of ref. 10. 
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the discrete spectra consisting of the harmonics of the modulation and its distortion (if 
any) due to over-modulation. When there is no signal, only terms in (3.6) and (3.7) for 
which m = 0 remain, and our general expressions of the present section reduce to the 
results of the preceding section, cf. Eq. (3.16) et. seq. Other limiting cases may be 
treated in the manner of section 4, reference [10]. 

Acknowledgment. The author wishes to express his sincere appreciation to Miss 
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ON THE NUMERICAL SOLUTION OF THE DIRICHLET PROBLEM FOR 
LAPLACE’S DIFFERENCE EQUATION* 


By 
J. B. DIAZ anp R, C. ROBERTS 
Institute for Fluid Dynamics and Applied Mathematics, University of Maryland 


1. Introduction. The present note is concerned with the proof of the convergence of 
three iterative methods for the numerical solution of the difference equation boundary 
value problem which is analogous to the classical Dirichlet problem for Laplace’s 
differential equation. These three iterative methods for the Dirichlet difference equation 
problem bear a marked resemblance to, and are patterned after, methods of solution 
for the Dirichlet differential equation problem, and may be briefly characterized as 
follows: (a) method I is the precise analogue of Poincaré’s [1]' “‘methode de balayage”’, 
(b) method II is the precise analogue of the extension of Poincaré’ method due to 
Kellogg [2], (c) method ITI is based on a certain connection between Laplace’s differential 
equation and the heat equation (intuitively put, ‘the steady state temperature is a 
harmonic function’’)’ 

Several convergence proofs for the Dirichlet difference boundary value problem have 
been based on Liebmann’s [5] original proof (see P. Frank and R. von Mises [4], and 
I. G. Petrowsky [5]) which involves a definite ordering of the points of the domain. It 
follows from the convergence proof given below for method III that this preliminary 
ordering of the points of the domain is unnecessary. R. Courant [6], has given a con- 
vergence proof for method I for the Dirichlet difference boundary value problem em- 
ploying a variational method which is the direct analogue of the classical Dirichlet’s 
principle for Dirichlet’s differential boundary value problem. In spite of the inherent 
interest and symmetry of the proofs in methods I and II described below, it should 
be remarked that the variational method of proof (minimization of a quadratic sum) 
proposed by Courant possesses a certain advantage in that it may be applied to other 
“elliptic” difference equation boundary value problems, since it does not require for 
its application the maximum-minimum property of the solutions of the difference equa- 
tion, i.e. that “the maximum and the minimum of a solution of the difference equation 
oecurs on the boundary’’. In the Laplace case, both the differential and the difference 
equations possess the maximum-minimum property. However, it is easy to give examples 
of elliptic differential equations (for which the maximum-minimum principle holds) 
e.g. U,, + U,, + U,, = 0, which have the property that the seemingly most naturally 
associated difference equation does not possess the maximum-minimum property. A 
convergence proof for method III, different from the one given here, has been given 
by Feller [7]. 

Although only the two dimensional case will be considered here explicitly, it is clear 
that the same arguments are valid in a finite number of dimensions. 

2. Convergence proofs. The boundary value problem will be formulated first. Con- 
sider a finite set S of points (m, n) in the zy-plane, where m and n are integers (i.e., S 
is “a net of equally space points in the plane, where for convenience the spacing has 


*Received May 4, 1951. This work was carried out under Contract N7onr-39705, sponsored by the 


Office of Naval Research. 
1Numbers in brackets refer to the bibliography at the end of the paper. 








356 J. B. DIAZ AND R. C. ROBERTS [Vol. IX, No. 4 


been taken as unity’’). The points of the set S are divided into two mutually exclusive 
subsets of points, the subset D consisting of interior points of S and the subset C of 
boundary points of S. A point (m, n) of S is said to be an interior point of S provided 
that its ‘‘four neighbors” (m + 1, n), (m,n + 1) also belong to S. A point (m, n) of S 
is said to be a boundary point of S provided that at least one of its four neighbors 
(m + 1, n), (m, n + 1) does not belong to S. It will be supposed in what follows that 
every point of the boundary C has at least one neighbor in D (this is clearly no restriction, 
as far as the difference equation problem is concerned, since only “isolated’’ points of 
the set S and “unessential’’ boundary points of S are excluded). 

The Dirichlet difference boundary value problem for the finite set S = D+ C 
consists in finding a real valued function u, defined on D + C, assuming prescribed 
values on the boundary C and which satisfies Laplace’s difference equation 


lu(m,n) = u(m + 1,n) + ulm — 1, n) + ulm, n + 1) + ulm, n — 1), (1) 


for each point (m, n) in the interior D. 

This boundary value problem, which is analogous to the classical Dirichlet problem, 
is well known to possess one and only one solution. This is readily seen as follows. The 
difference equation (1) amounts to a system of d linear equations in the d unknown 
values of u at the points of D, supposed to be d in number. This system of linear equations 
is homogeneous if and only if the prescribed values of u on C are all zero (recall that 
it was assumed that every point of C has at least one neighbor in D). But it is known 
from algebra that a system of d linear equations in d unknowns will have a unique 
solution if and only if the corresponding homogeneous system has only the trivial, 
identically zero solution. That the homogeneous system has only the zero solution 
follows immediately once it is shown that the maximum and minimum values of a 
function u, defined on D + C, and satisfying (1) in D, must occur on the boundary C. 
Suppose, on the contrary, that the maximum value M of u occurs at a point (m, n) of D, 
but not on the boundary C. Then from (1) it follows that the value of wu at the four points 
(m + 1, n), (m, n + 1) must also be M and hence these four points must all be in the 
interior, D. By induction, it follows that the infinite set of points (m + j, n), 7 = 0, 1, 2, 

all belong to D, which contradicts the initial assumption that D is a finite set of 
points. Thus the maximum value W/ must occur at a point of C. Since a minimum of 
the function uw is a maximum of the function —u, and —w also satisfies (1), it follows 


that the minimum of uw on D + C also occurs on C. 


Method I. 

Consider the convergence of method I mentioned in the introduction. Denote by 
f the given real valued function, defined on C, which is the prescribed boundary value 
of wu on C, and let the “initial function” G be a superharmonie function defined on D + C 
and which coincides with the given function f on the boundary C. That is 


bo 


1G(m, n) 2 Gm + 1, n) + G(m — 1, n) + G(m, n + 1) + G(m, n — 1), ( 
for each (m, n) in D, and 
G(m, n) = f(m, n), 
for each (m, n) in C. The restriction that the initial function G be superharmonic is not 
essential, and will be removed later. 
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Let the points of D be arranged in a sequence P, , P, , P; , --+ in such a way that 
each point of D occurs infinitely many times in the sequence. For convenience, the 
sequence may be chosen by moving from point to point of D in some definite geometrical 
pattern, but this is not necessary. 
A sequence of functions is now defined on D + C by the following procedure. First, 
let wy G, and define the sequence of functions wy , w, , W, , «++ successively as follows: 


First step: w = G, on D + C; 


Second step: w, is harmonic at P, , 


w, =wonD+C—P,; 


Third step: wy, is harmonic at P, , 


w, =w,onDdD+C-P.,; 


l 


p + 1) st. step: w, is harmonic at P, , 


w, = w,,onD+C-—-P 


yp? 


etc. In other words, one moves from point to point in D changing the value of the func- 
tions at each point so that equation (1) holds. 

It is easily seen that each function w, is superharmonic in D + C. Consider w, , 
and let P,; = (m, n). w, is certainly superharmonic at all points of D + C, save perhaps 
at the five points (m, n), (m + 1, n), (m, n + 1). But at (m, n) the equality sign holds 
in equation (2), with G replaced by w, , by the definition of w,(m, n), while at the four 
points (m + 1, n), (m, m + 1), equation (2), with G replaced by w, , holds a fortiori, 
since the right hand sum does not increase if G is replaced by w, . By induction, w, 
is superharmonic for each p. Further, the sequence of functions wy , w, , w. , «+ is 
monotonically non-increasing on D + C, by the way in which the sequence was con- 
structed. Finally, each function w, is bounded below by the minimum value of G on 
D+ C. 

This last statement follows immediately from the fact that the minimum value of a 
function superharmonic in D + C must occur at a point of C (this may be seen by an 
argument analogous to that used above in proving the maximum-minimum property 
of solutions of (1)). Since 

w,(m, n) = G(m, n) = f(m, n), 
for (m, n) in C and any p, it follows that for any p 
w,(m,n) 2 ming G = ming f, 
for any (m, n) in D; and this, together with 
w,(m, n) = w.(m,n) 2 --+ 2 w,(m,n) 2 --- 


assures the convergence of the sequence w, at each point of D. It only remains to show 
that the limit function 
lim w,(m, n), 


p72 


(m, n) in D + C, is the solution of the boundary value problem. 
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There is no question about the fact that the limit function coincides with f on C, 
since each w, does, it is only necessary to show that the limit function satisfies the partial 
difference equation (1) in D. Consider a point (m, n) of D. Since (m, n) occurs infinitely 
many times in the sequence P, , P, , P; , --+ it follows that there is an infinite sequence 


of integers, k, , such that 


(m,n) = P,,, 
for p = 1, 2, 3, --- . Consequently, there is a subsequence of functions w,, , w,, , ; 
Wr, , °** Of the sequence w, , w., +++, W,, *:+ Such that each function w,, of the subse- 


quence is harmonic at (m, n), that is 
4w,,(m, n) = w,,(m + 1, n) + u,,(m — 1, n) + wy, (m, n + 1) + w,,(m, n — 1). 
Since 


lim w,(x, y) = lim w,,(2, y), 


M poo 


for any (x, y) in D, it follows that 


4 lim w,(m, n) = lim w,(m + 1, n) 


po x 


+ lim w,(m — 1, n) + lim w,(m, n + 1) + lim w,(m, n — 1), 
ou ie oe 
i.e., the limit function is a solution of the boundary value problem (notice that the 
existence of a solution has been proved independently of the earlier considerations in- 
volving Cramer’s rule). The solution has already been shown to be unique by the 
maximum-minimum principle. 

The restriction that the initial function G be superharmonic will now be removed 
by showing that any function defined on D + C may always be represented as the differ- 
ence of two superharmonic functions on D + C. Let g be any function defined on D + C, 
and let .V be one fourth of the maximum of the absolute value of the Laplacian of g on D, 


i.e. WM = 4+ maxp | g(m + 1,n) + g(m — 1,n) + g(m,n+ 1) + g(m, n — 1) — 4g(m, n) 
Then 

g(m, n) = gi(m, n) + go(m, n), 
where 


gi(m, n) = g(m, n) — M(m* +n’), 


g.(m,n) = —M(m’* + n’). 


Clearly, g, and g, are both superharmonic on D + C. Thus, if g is any function taking 
the prescribed values on C, the sequences of functions, constructed as explained above, 
corresponding to the initial superharmonic functions g, and g, converge to harmonic 
functions, and the difference of the two limit functions is the desired solution. 

The method just described is the precise analogue of Poincaré’s [1] “method of 
sweeping out” for the classical Dirichlet problem. It also coincides with the relaxation 
method (Southwell [8]) provided that the relaxation is done point by point and that at 
each point the residual is actually reduced to zero. 
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Method II. 

This method differs from the preceding one only in the removal of the restriction 
that the process of method I be carried out pointwise at each step. Consider a sequence 
B, , B, , B, , --+ of subregions (“blocks”) of D, subject to the following conditions: 
(a) each point of D is an interior point of an infinite number of the blocks of the sequence, 
(b) the Dirichlet problem is explicitly solvable for each block B, in terms of arbitrary 
boundary values on the boundary of B, . 

The convergence proof is the same as in method I. One need only substitute the 
sequence B, , B, , B, , --+ of blocks for the sequence P, , P; , Ps , --+ of points. 

This method is the exact analogue of the procedure suggested for the classical 
Dirichlet problem by Kellogg [2]. It is not to be confused with “block relaxation’, with 
which it coincides only in the case where the residuals are actually made zero at each 


interior point of each block. 


Method III. 
Method III is essentially a method of successive approximations. Starting with an 


initial funetion w) = g which satisfies the prescribed boundary condition, one defines 


the following sequence of functions w, : 
Ww = g,onD+C; 
4w,(m, n) = w(m + 1, n) + w(m — 1, n) + wom, n + 1) 
+ w(m,n — 1), for (m, n) in D, 


w,(m, n) = g(m, n), for (m, n) in C; 


tw,(m, n) = w,-.(m + 1, n) + w,-.(m — 1, n) + w,-1(m, n + 1) 
+ w,-i(m, n — 1), for (m, n) in D, 


w,(m, n) = g(m, n), for (m, n) in C; 


Again one may suppose, without loss of generality, that the initial function g is super- 
harmonic to start out with. Hence all the functions w, are superharmonic and from the 
minimum principle for superharmonic functions it follows that 


Wy SW, > WwW. See? Sw, = ees JS Ming Gg = Ming f. 
i= 2 Dp e d 


Thus the sequence converges, and since 


lim w, = lim w,-; , 


p72 pa 
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it follows that the limit function of the sequence is the desired solution of the Dirichlet 
problem. 

The last method is related to the difference equation formulation of the heat flow 
problem, Emmons [9], in which the boundary C is held at all times at a fixed temperature, 
g represents the initial temperature, and the functions w, represent the temperature of 
D + C at regular intervals of time after ¢ = 0. The limit function represents the steady 
state temperature. A proof of the convergence of method III, based upon random walk 


considerations, has been given by Feller [7]. 
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ON THE NON-LINEAR VIBRATION OF ELASTIC BARS* 


BY 
A. CEMAL ERINGEN 
Illinois Institute of Technology 


1. Introduction. The classical theory of vibration of bars is based on certain restrictive 
assumptions, namely: (a) the deflection is small; (b) supports are free to move in the 
axial direction; (c) deflection is inextensional. In practice, however, very often some or 
all of these assumptions are violated. Therefore, it is necessary to reformulate the 
problem of vibration of bars in its general form without these assumptions so that the 
domain of applicability of the classical theory can be well defined and problems to which 
the classical theory is not applicable can be attacked. 
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Fic. 1. Displaced element of beam. 


Recently Woinowsky-Krieger [1] studied the effect of axial force on the vibration of 
hinged bars. N. J. Hoff [2] gave an analysis for the effect of inertia forces on the buckling 
of columns. In both of these analyses the classical treatment was improved by considera- 
tion of the axial stress due to bending. 

In the following analysis, first the problem of vibration of bars is reformulated 
without any of the above-mentioned assumptions. Following that, a basic problem, free 
vibration of elastic bars having immovable hinged ends, is solved with the use of the 
perturbation method. The solution of this problem adequately describes those motions 
in which the changes in axial tension, as well as in deflection, are not small. 

Solutions of vibration problems concerning bars with other types of end conditions 
and foreed vibrations are left to a forthcoming paper. 


*Received Feb. 15, 1951. 
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2. Equations of motion. The equations of dynamic equilibrium of an element of the 
beam deformed to a plane figure are (Fig. 1) 


—pAu.., + (N cos 6), — (Qsin @),, = 0, 


I| 
S 


—pAv...+ (Nsin 6), + (Q cos 6),, 
(1) 


(1 +u,.) J0..+ M,. — Q[1 + u,,) cos 6+ 0,,8Iin 6] 


+ Niv,, cos 6 — (1 + u,,) sin 6] = 0. 


where p is the mass per unit volume, A the cross-sectional area of the beam, J the mass 
moment inertia per unit length, NV, Q, M are the axial force, the shear force, and the 
bending moment at a point z, u and v are the deflections of any point z in the axial and 
transverse directions, @ is the angle between the tangent to the median line and the 
x axis, and / is the time. In calculating the inertia forces, the effect of changes of mass 
during motion is neglected. 

The length element, ds, after deformation is given by 


ds’ = {[v.. + (y cos @),.) + [1 + u,. — (ysin 6),.]'} dx’, tan@d=v,/(1+u,,). (2) 


— t 
Indices after a comma represent differentiation. In deriving Eq. (2) the effect of shear 
deformation is neglected. 

Normal strain is defined as 


ds 
=—-— |] (3) 


é —_ 
dx 


Equations (2) and (3) are combined to yield 


BS Sa 2 op St.» I, (4) 


e=e— y8, é= hari 
YP.z cos 6 sin 0 
Here, « is the strain referred to the median line of the undeformed beam. 


Hooke’s Law states that 


Thus, 


’ 


N=|[cdA=EA, M= | oy dA = —EI@.,, (6) 


are calculated with the use of the first Eq. (4). Here J is the moment inertia of any 
section about its neutral axis. The variation of J with time is neglected. 
Equations (6) and the third Eq. (1) are combined to give 


Q= ———S + J6 4, cos 6. (7) 
The first two Eqs. (1) can be transformed into more suitable form by introducing 


non-dimensional quantities and multiplying the second by 7 = ~/—1 and adding the 
result to the first. The resulting equation is differentiated with respect to x. Thus, 
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—[(1 + ee]... + {[. + a(— ee + 6.,, cos a) ket = 0, 
(8) 
r = (E/pL’)'2, y = 2/L, ? = J/pAL’. 
The second Eq. (4) may be transformed to 
y + (u(x, t) + iv(a, t) — uO, t) — w(O0, d] = | (1 + ee" dy. (9) 


The first bracket in the first Eq. (8) represents the effect of translational inertia; 
the first, second, and third terms in the brace are respectively the contributions of the 
extension of the median line, the bending, and the rotatory inertia. 

For \ = 0 the first Eq. (8) reduces to the one obtained by Carrier [4] for the problem 
of non-linear vibration of the elastic string, when ¢ is replaced by (€ — €) in order to 
include the initial tension. 

Further, let 


s = dr = U(EIg/WL*)"”, - 


where IV’ is the total weight of the beam. Equations (8), when separated into real and 


imaginary parts, become 


2n°6 [ 2 €.,0.,0 v 











—n*e.. + (1 + €)d76", ,— 68, - So — 2? Sw 
€ + ¢ + €y «0, + er r (1 + 
(11) 
+X a 2r°6...,8,, cos 0 + 2r°O,,,0, sin 6 — d°O,,,8,,, cos 6 = 0, 
- € 
; . O.. 2 €,0 
—2r'e. 0, — (1 + ors, I~7- i 2° eet, 
€ € es t 2€,0., + €6.,, r 1 7 (i + 6? 
—_ : ee >» 6°,0 ‘ ‘ 
+ *- ww — 2°—*45 ~ — 0 8Y + "8 seyy COS 8 — 2A'O.,0..2,8mM B 12 
(1 + 6&) ica’ ite* cory COS sis (12) 
— 2\*6,0.,, cos 0 — d*0,,,0.,, Sin @ = 0. 
Contrary to expectation, it can be seen that « = 0 is not a possible solution since 


Eqs. (11) and (12) reduce to two independent equations for @. Consequently, assumption 
c) of the classical theory is not valid. It can be seen, however, that this assumption 
will be valid when « is a higher-order quantity than 6. Therefore, the classical theory is 
obtained as a limiting case of the present theory for @ — 0. 

3. Beams with immovable hinged ends. The boundary conditions for simply- 
supported beams with immovable hinged ends are 


6, = 0, u=v=0 for y = 0, 1. (13) 

The perturbation parameter is chosen as \ which represents the ratio of rotational 

inertia to translational inertia. Reversal of the sign of \ should reverse the sign of 6 
but not e; @ and ¢ are therefore expanded into power series of \ as follows 


6= 0, + r\°6, + XO, + ---, e=Ne+r‘at+°::. (14) 
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The expressions for @ and ¢ are substituted into Eqs. (11) and (12) to obtain differential 


equations and into (13) and (9) to obtain boundary conditions. 


Differential equations 





(X”) F = 0) 
(r*) €4 = € = 6; 20 7] o—_ 6; y = %e € + €,0; vy?) 
(X” 15) 
(r° 6 + 6 Oe, = Bath. « =O, 
(n°) 6 + 9 — €65 = 6 — 2, .0,., — 460, 
2.8; + 2e, 6; + ¢, ,,6 + 6; 0; + 2, ,6 
Ze, ,0,;., + Geves ,6;., + €49; + 356, = ( 
(X‘): ane % 16 
Bounda 1 condition 
at y 0, 1 6 = () 1=1,2 ) li 
+) A 
Ae OF Cede 0, nN | 2 ) dy = §, 
, A 
» ] - A if 4 dy 0, | € P ia 6,6. = 2 - ) dy ( 
\ 6 \ 9 9 
Ss : Sy, (A jt =** 3 19 
Tnitial condition 
at S (): i] mré,, COS mMrYy, A 0), } Q € 0 
() 
Fee ee 
The solution of the first hq 15) under the condition that €.(0, s) €,(1 and 
the first Eq. (19) is 
. ) 
€ | 0 dy ] 
Substitution of (21) into the first eq. (16 leads to 
”») 
A + 6 6 | 6; dy Q) 22 
2 P 
Equation (22) is solved under the boundary conditions (17) and the first I: 8 
for an initial sinusoidal deflection given by the first Inq. (20). 
Let 
mrG@,- COS Mri, S(O l. ae 
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gives the following differential equation: 


Substitution of lq. (23) into Iq. (22) 

+ (mar)'S + 1(mr)'6S* = 0. (24) 
Multiplving (24) by S a first integral can be obtained immediately. The resulting 
equation is of first order and separable. Integration can be effected in terms of elliptic 
functions. The result is 


7 fon 
0 MmMrg, COS MmMrY CN (WS, k), oO, = m (1 + \ 4, ; 


; \ "yy 9 2\-1/2 
p= (245) , Ps = 2(1 + 8) K, 


vhere cn and A denote the elliptie cosine and the complete elliptic integral of the first 
kind; 7, 7 


is the ratio of non-linear period to linear period. The use of Eq. (21) gives 


(26) 


€, = 1(mr6,)” en” (w, 8, k). 


Substitution of (25) and (26) into the second Eq. (15), integration, and the condition 


€.(0 é. 1, s), vield 
€ e(s) cos 2mry + e,(s), 
4 
mn) |g? +! 6: + 26° en’ | k) + St wal Gece: (27 
8 lL — @, 26, en” (ws, k - #, en (w,8, &) |, 27) 
s [ 8 te 8 
5 +1 
» ‘ 4 P 
(s m6.) en’ (ws, k) + mA, en (ws, /) | 6, cos mary dy, 
{) : “0 
where is obtained after the use of the second Eq. (19). The second Kq. (16) thus 
he« les 
I ° 12 2 
i] 1 # ee mar @ cn” (ws, h)@ + (mm) 6, en” (ws, k). 
CO Tu |) Os. COS mary dy Rj(s) cos mry + R.fs) cos 38mr7y, 
y 12 13 15 
| ” 2) ” ; 
he mr)'| (1 —-— 6 0) 6, en (ws, k) + = 65 en” (ws, k) 
q 16 | ae Ss 
(28 
21 
+ 4, ¢ (ws. K) 
i 128 cn @,8, / Be 
) 5 3 
R yar) | 6 en (ws, k) — 6; en (ws, k) — = 6, en” (as, k) 
It) | ) g 
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The solution of partial differential equation (28) is of the following form: 


6, = 0,” + 0," + 43°, (29) 


where @;” is the solution of the reduced equation. 6@;'’ and 6;°’ are particular solutions 
of the differential equations obtained by taking F,(s) 0 and R,(s) = O respectively. 
An examination of the second Eq. (18) reveals that 


6,” = S,(c) cos mry, 6, = S.(c) cosrry, og = oS, r~ m, (30) 


where S, and S, satisfy the following differential equations which are obtained by 
substituting Eqs. (30) into the reduced equation 
g 


l,ooe 


= [2k* sn? (c, k) — 1]S,, (31) 


So.ce = [6k’sn* (co, k) — (1 + 4k°)JS, . (32) 


Differential equations (31) and (32) are the Jacobian forms of the generalized Lamé 


equation. The solutions of these equations are: 


S. = BA,’ +CA.’, (a = 1, 2) (33) 
AS’ cen (a, k), E 
(34) 
A; = en (oe, k)[U1 — k)o + dn (oe, k) se (0, k) — Eo, k)); 
= H(oe + a,) P 
a : “~ exp [—oZ(c,)], 
; Sere p | 
= 
(30) 
; + H(c — a, a 
A, = I] ——* } exp [oZ(c,)], 
; O(c) 
where ¢, and a, are chosen to satisfy the following two independent equations: 
sn og, cn o, dn o, + sn oa, cn o2 dn og = O, 
36) 
en o, ds o, + en a, ds os)” — ns o, — NS a» —(1 + 4k’). 


Here @(u), H(u), and Z(u) are Jacobian Theta, Eta and Zeta functions respectively [4]; 


E(u) is the fundamental elliptic integral of the second kind; sn u, en u, and dn u are 
Jacobian elliptic functions and sc u = sn u/en u, ds u = dn u/sn u, ns u = 1 sn was 
originated by Glaisher. These functions are tabulated in [6]. 

The general solution of an equation similar to (32) was first given by Hermite [5]. 
Only the first of Eqs. (34) can be extracted from this solution by use of properties of 
these functions. The second solution A;” is obtained upon reducing the order of the 
differential Eq. (31) by letting A,” = A(c) en (o, k) and upon integrating the resulting 
equation for A(c). It can be seen that while A;’’ is a doubly-periodic function A,” is 
not periodic. 

Particular solutions 63’? and @;°’ of Eqs. (28) are of the following forms: 


6," = S"’(e) cos mry, 6; = S*’(c) cos 3mry. (37) 
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Substitution of Eqs. (37) into the first of Eqs. (28) with R,(s) = 0 and R,(s) = 0 


respectively, results in ‘ 
S% — [2k’? sn’ (c, k) — 1]S = R,(o), (38) 
S°) — [6k? sn? (o, k) — (1 + 4k*)]S™ = R,(o). (39) 


Particular solutions of these equations are obtained by Lagrange’s method of varia- 
tion of parameters. Hence* 
1 ag 


So) = 5 | (AL) APH — AP HALON at, (40) 


4=-(1-F), @=1,2) (41) 


> 


A. = 2=— ye — cn o; dno; nso; 42 
2 6°(0) > o; dno; nso; , (42) 
where o, and o, must be solved from Eq. (36). Finally, the use of the initial conditions 
given by Eq. (20), namely, @; = 6;,, = 0 for s = 0 reduce the general solution @;(y, o) to: 

6,(y, 0) = S‘(¢) cos mry + S'(c) cos 3mry. (43) 


e,(s) of Eq. (27) thus becomes 
3 ae 1 v1) 
e,(s) = GA (m6) cn’ (ws, k) + 9 mr OoS-°'(w,s) en (ws, k) (44) 


This completes the solution of the equations corresponding to \* and d°. Further approxi- 
mation requires more tedious analysis. However, since larger initial deflection will cause 
extremely high axial tension in the bar, the validity of Hooke’s Law must be examined 
before taking any further steps toward improvement. 

There still remains the examination of the question of convergence. This is impossible 
at this point. Although a term, ¢, appears in the expression for A,~’, this does not neces- 
sarily mean a resonance effect, since further approximations are necessary to examine 
the series in o for the question of stability. 

An estimate of region of stability of Eqs. (31) and (32) can be made for small &. 
In this case sin o can be used in place of sn o. Therefore Eqs. (31) and (32) become 
Mathieu equations whose stability regions in terms of parameter k* are well known [7], 
and can easily be seen to contain those of the present equations. 

4. Vibration produced by an arbitrary initial deflection. Analysis of the preceding 
section is based on an initial condition corresponding to an initial sinuzoidal deflection 
of the hinged bar. Differential equations must be re-solved when the initial deflection 
is arbitrarily prescribed, since these equations are non-linear. The problem of non-linear 
vibrations of the elastic string following an initial deflection has its analogue in the case 
of hinged elastic bars. Following a similar method, as that of [3] an integral equation is 
obtained below which may be solved by the method of successive approximations. This 


*Eq. (40) for a = 1, can be integrated in closed form, in terms of elliptic functions. This result will 


not be given here, however. 
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equation represents the A°-approximation of the problem since it is obtained by integrating 


or 22 Let 
ing 
A(y. 8 - | (8S) COS mMTY, 
(45 
A(y4. O pe | COS MT, S,,(0 i 


where | are the Fourie coefficients for the initial deflection function @ i, 0 Sub- 
stitution of Eq 15) into hq 22) leads to 


' (mam) — ‘ 
S, T (mr) SS, ~ s > As. (46 


The required result is obtained by integrating Kg H)): 


; Fm bes 5 
S,(8) = COS mrs — ’ | sin mam(s — O)S,(0 Zz. AiS,(t) dt. (47 
There remains the problem of solving €, , &; , «++ : 6: . 6. . et , for an arbitrary initial 


deflection. However, corresponding differential equations given by (15) and (16) are 


linear. Therefore superposition is valid after 6,(y, s) is solved from non-linear integral 
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Fig. 3. Axial stress versus period. 


It is seen that axial stress increases very rapidly with the decrease of 7/7) or with the 


Increase ol frequency ratios. 
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ENERGY THEOREMS AND CRITICAL LOAD APPROXIMATIONS 
IN THE GENERAL THEORY OF ELASTIC STABILITY* 


BY 
J. N. GOODIER anp H. J. PLASS 
Stanford University 


1. Introduction. When the ordinary uniform pinned-end column buckles under a 
critical load P, = x°EI/l’, the potential energy measured from the straight compressed 
form is 


ae 1] ae 
V =-EI '? dy — =P a | 
or 
. al 2 I 
io =| yde—%] y? de (1) 
4 “0 “0 


and is zero since Y « sin wa/l. 
Let y now be a deflection of any form which satisfies the end conditions. Then V 
\ 
Po k—-_ ¢ ——_+ P 





Fic. 1. 


as given by (1) is positive unless y « sin ra/l, in which case it is zero. This follows readily 
from a Fourier expansion’ of y, or from Wirtinger’s inequality,” which is that 


a25 ale 


| 2 dx < | 2’? dz, (2) 
unless z = A cos x + B sin 2, provided that 
2(0) = 2(27) and [- zdz = 0. 
Jo 
For if we reflect the bar with its deflection y in its right hand end (x = 1) and combine 


the inverted reflection with the original bar (Fig. 1), we have a bar of length 2!, with 
y'(0) = y'(2l), and from (1) 


4V an r72 a ss 72 
EI > yy’ dz — E y’” dx. 


*Received February 5, 1951. 

See for instance R. V. Southwell, An introduction to the theory of elasticity for engineers and physicists, 
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The substitution — = ra//l converts this to 


“ 


where the prime now means differentiation with respect to & and if we identify y’ with 
z, € with x in (2), (2) establishes that V is positive unless y is sinusoidal. 

Two conclusions follow. Firs/, that since V is the potential energy of the bar under 
the true critical load P, , the sinusoidal buckled form is itself stable with respect to 
disturbances (impulses) which project it into non-sinusoidal forms. It is neutral with 
respect to sinusoidal disturbances. In the disturbances ?,; remains unchanged. Second, 
that an approximation P, calculated by inserting an assumed approximate deflection 


vy, tor yin the relation 


arr aor 
9 EI yo dx — sf y Yo ax 0 (4 
valid hen P? is P nd y is sinusoidal) will be higher than ?, . For the use ot 1) in 
this v will vield 
EI | yi? d 
P 5) 
| y, d 
But sinc is not sinusoidal we have > and (1) vields 
EI | yf? 
r 6 
| y ad 
Here e have ! nequality 2) available from pure mathematics, and ean use it to 
establist either tne first o1 the second conclusion 
iF plate problem the corresponding inequality, establishing that the potential 
energ\ strain energy ol bending minus work of critical loads on buckling displace 
{ any displacement differing from the true buckling displacement 


ments), is positive for 
is a much more elaborate one, although it can of course be formulated. A proort DV means 


of Fourier series feasible for the simple) cases, such as the rectangular plate vith fou 
simply-supported edges, but a proof for the more difficult cases such as four clamped 
edges is h irdly to pe expected Still less can we hope to obtain such a proort tor more 
complex systems such as shells, or combinations of structural elements such as stiffened 
plates and shells, or the general problem of elastic stability with respect to infinitesimal 


displace ments 
Sut if we are given that the buckled form is itself not unstable, this datum establishes 
the inequality, and we can then use it to prove that the energy approximations to the 
critical loads will be too high. In the remainder of the paper we do this for the general 
stability problem. If the buckled state were itself unstable, the energy approximation 
to the critical load would be too low. Thus the usual assumption in practical calculations 
that the approximation will be too high is equivalent to the assumption that in the 


idealized version of the problem there 7s a buckled state which is itself not unstable. 
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2. Formulation of the general equations. An arbitrary elastic solid has initial stress 
specified by the usual Cartesian components S,; (7, 7 = 1, 2, 3), which maintain equili- 
brium with initial body foree F, per unit volume and surface force 7’; per unit area on a 
surface element whose outward direction cosines are v,. We have then the differential 
equations of equilibrium (with the summation convention for repeated indices, and sub- 
scripts after a comma indicating differentiation with respect to the corresponding co- 


ordinates 


S,, + F, =0 7) 

and the boundary conditions of equilibrium 
Sv, = 7 S) 
The stress S,, is not necessarily entirely due to F, and 7, . It may be initial or thermal 


tress existing In the absence of fF, and 7’ 
‘ss Will be referred to as state 1, and x, are the co-ordinates of material 


Chis state of stre 
points in this state (not in the unstressed state). For the present, we suppose that it is 
stable. A second state, state II, is derived from it by the application of additional body 
force AF, and additional surface force AT, . The displacement caused is expressed by 
Cartesian components u, (not Au;), and it is affected by the presence of the initial stress. 
Phe stress in state IT is of course different from S,; . To specify it we use Trefftz’s stress 


components” / in Itappus’ notation). These are non-orthogonal. A rectangular block 


element in state [ becomes an elementary parallepiped in state II, and these stress 


components refer to the directions of its edges. The advantage of using them is that they 


lead to relatively simple equations. We may write 

/ S e- ¥ 9) 
und 7 r,, since both S S,,andk,; = k,, . Even where the 7,,; vanish, this state 
of stress need not be identical with that expressed by the S,, of state I, on account of 


the different specification of stress components. 


The differential equations of equilibrium’ satisfied by 7,, are 


riii + (Sati), J + OF; = 0 10) 
after neglecting “non-linear” terms (7,,u;.,).; , and so restricting the investigation to 
small compared with S more precisely to the largest 7,, small compared with the 


largest SS 
The boundary conditions of equilibrium satisfied by 7,, are 


‘i + 8H:.@; = AT; . 11) 


When S 0, (10) and (11) reduce as they should to the equations of the ordinary 


theory of elasticity 


equations (10) and (11) do not involve any stress-strain relations. Being concerned 


with small departures from state I, we assume that small strain components 
e,; = Huy; + u;.;) 12) 
Trefitz, Zur Theorie der Stabilitdt des elastischen Gleichgewichts, Z. angew. Math. Mech. 12, 160 


See the reference in footnote (3). 
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are related to the small stress components 7,;; by the usual form of Hooke’s Law. It is 
convenient to take this in the general form appropriate to the homogeneous anisotropic 
solid, and to express this form in a changed notation. Write 7,,; = 7, , T22 = T2, 733 = Ts, 
Ti2 = Ts, T23 = Ts, T3, = 7. and similarly for the strain components. Then the stress- 


strain relations are” 
tT; = ¢;;e; , Where 7 and j have the range 1 to 6, and ¢;; = ¢;; (13) 


The variational principles (stationary potential energy for equilibrium, Castigliano’s 
Theorem, ete.) of the ordinary theory of elasticity can be derived by considering the 
variation of the strain energy of the body, and using the equations of equilibrium.® We 
follow this method now for the transition from state I to state II, with some modifica- 
tion, using the equations of this article. We thus regard equations of equilibrium (or 
motion) as basic, and energy principles such as stationary potential energy as derived, 
rather than vice-versa. 

3. A variational principle. From the quantities e,(¢ = 1 to 6) as functions of the 
co-ordinates z;(7 = 1 to 3) in state I we may form by integration over all the volume 
elements dw of state I the integral 


ny 


U(e) = >| C;;ee; dw (14) 


which would be the strain energy in the absence of initial stress. Trefftz’ has shown that 
the'strain energy acquired in the passage from state I to state II is 


a 


: ss is . ne : ss 
U(e) + | S;e;; dw + = | S;,U;, Ui.- do, (zt, 7 = 1, 2, 3). (15) 
We consider the variation of U(e) alone. Let arbitrary variations é6u; be added to the 
displacements u; . Then, writing 6U for the complete variation of U we have from (14) 
= — : + i 2 bos : 
U(e) + 6U = z c;;(e; + 6e;)(e; + 6e;) (7,7 = 1 to 6) 


a? 


YF < mae, 
= U(e) + | c;;(e;6e; + e,;6e;) dw + U(ée), 
where 


- EF . 
U(ée) = | c;; 6e; 6e; dw. (16) 


Since ¢;; = ¢;; , we have 


| 7; 6e; dw + U(ée). (17) 


d 


6U = [ cisesde, dw + U(ée) = 


We now return to the range 1, 2, 3 for 7, 7 and k, and write instead of (17): 


6U = | ris 60:5 dw + U(6e). (18) 


‘ 


5A. E. H. Love, Mathematical theory of elasticity, Cambridge University Press, 1927, Ch. 3. 
6. Trefftz, Handbuch der Physik, Vol. 6, Springer, Berlin, 1928, p. 68. 


7See the reference in footnote (3). 


o 
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For the integral in (18) we have, taking ée;; from (12), 


Ls : 
: | 7; ;(6u;.; + 6u;.;) dw = F% 6u;,; dw (since 7;; = 7;;) 


= [ 6u;),; dw — Cc bu; dw. 


By an application of the divergence theorem to the first of these two integrals we may 


write the result as 


7;; 6u;v; do — / 7:i;,; 6u; dw, 
where do and & refer to the boundary surface and v; are the direction cosines of the 
normal, all in state I. We now eliminate 7;; by use of the boundary conditions (11) in 
the first integral and the equilibrium equations (10) in the second. The result is 


[ [—Sjusar; + AT.) bu, do + | ((Sjatts.a).; + AF ;] bu, deo (19) 


The first part of the second integral is transformed as follows 


| (S;,u;..).; du; dw = [ (Sins. bu;),; dw — il S5.U;,; bu;,, dw 


oe [ S;uU;.. 6u;v; do — [ Sian. du, dw. 


With this (19) is simplified by cancellation of the two surface integrals involving S;, . 
Recalling that (19) is equivalent to the integral in (18), we have as a new version of 


the latter equat ion 


6U = | AT; bu; do + | AF; bu; dw — [ S;,u;,; 6u;,, dw + U(ée). (20) 
In this the u; are the actual displacements caused by the application of the additional 
forces AT, and AF; , corresponding to the passage from state I to state II, and the 
du, are arbitrary additional displacements. Both the u; and the 6u; are restricted to 
smallness by (12) and (13). 

Now let S,, , F; , JT; , AF; , AT; be fixed, but let wu; for the moment be three inde- 
pendent functions of the x; , not required as yet to be the correct displacements in the 
passage from state I to state II. The result (20) suggests consideration of a function of 
these u,; in the form 


. , 


, es oa J 1 = 
V = U@) — | AT ju; do — | AF ju; do + > | S51; Ui.n dw. (21) 
On varying the u; (as they appear explicitly, and also in U(e)) we have 


6V = 6U — | AT; du; do — [ AF; 6u; dw + ; | S;, 5(u;, U;,,) dw. (22) 


Jy 
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Now let w. be eiven the actual values of state IT. Then on account of (20) we have 


; V [ S : 6(u d . { bu dw -t- U( Or ). 23) 
sult 
iS i / cS uw. + Ou.) (Ul + i — ul 
aS | Ou u;4 + u Ou 1 6 ou ¥ 
Since S S he first term in the brackets can be combined with the second to give 


the result 


LS ..12u Ou 1 Ou bu 
Then (23) reduces to 
. Sf ae 

5 | QO+ U6 - | Ss Ol Ou dw. 24 
the zero indicating that the first order (in é6u variation of Vo vanishes. This property 
of course would be characteristic of the potential energy in state I] as an equilibrium 
state. Rete rred to state I as zero, the potential energ\ consists of the strain energy 15 
togethe vith the potential energy Ol the body and surface forces, which is given for 


It can be shown by means of (7) and (8) that the terms here in 7’, and F,, cancel the 
middle term in (15), and hence that V as given by (21 is in fact the potential energy 
of state IIT when the w, in (21) denote the actual displacements of state IT. 

4. The stability of state II. Our object being to deduce a generalization of the in- 
equality (6 when the stability of state IT 1s given, we now seek a necessary condition 
for this stability 

Let the particles of the body be projected from state [1 by some disturbance. Then 
at time / they are in motion with displacements éu, (functions of the co-ordinates x, of 
the particles in state I), and corresponding to 6u; we have additions 67,; to the Trefftz 
stress components 9). There is, by hypothesis, no change in the forces fF, + AF 
T, + AT, of state II except at fixed supports where reactions may be induced. During 
the motion these are carried with the particles on which they act in state IT. The equa- 


tions of motion are, from (10) 
- + 67...) aa 1S: u : bu.) xl - AF = pou, ’ 


where 6ii, = 0°6u;/dt (the acceleration) and p is the density. Subtracting (10) we have 


Tii.i + (S,,.du DY i pou 


Multiplying by 6u,; and integrating over the volume we find 
F Fa Pe 3 d 
T 6u; dw + CS; 6u ) ; Ou; dew = 


© 


Bat al 
7 | 5 P OU; Ou; dw, 2) 


1952 GENERAL THEORY OF ELASTIC STABILITY 377 


and the term on the right is the time derivative of the kinetic energy. The first integral 
on the left of (25) can be written as 


> 


| (67;; 6u;).; dw — | 67;; 6u;.; dw 


“ 


and after transformation of the first of these integrals by the divergence theorem, as 


> 


| 67r,; 6u,; v; do — 67;; 6u;.; dw. 


. 


Sim lv the second integral on the left of (25) transforms into 
| S.. 6u 6u,v, do — | S.. du, , du dw. 
Introducing these transformations in (25), and writing 7 for the kinetic energy, we have, 


le rearrangement, 
I (6r,; + Sy Su;.4) 6u;,; dw = | (67,4; + Six bus) bu, v; do. (26) 
he cket appearing in the integral on the right is, by (11), the addition to AT; ac- 


ne the motion, and this is zero by hypothesis except at fixed supports, where 


ish. The integral therefore vanishes. The integral on the left of (26) is the 


 - i 7 : _ 
| OT Or dw + | Sop OU; Ou dw ?. (27) 
dt \2 i 
Th = readily erified by changing 67, ; 6¢ to 67,6e,(¢ = 1 to 6), then to c,;6e,6¢ 
carrving out the differentiation with respect to ¢, and making the combinations of terms 
permitted by « c,, and S S,. . It is evident that the bracket in (27) is identical 
las given by (24). We can therefore re-express (26) as 
IT’ 16V 
a ado 
ee i - = 0 (28) 
at dt 
howing that 7’ + 6) remains constant during the motion following the projection from 
the equilibrium state II. This of course is the energy equation of this motion, and 


exhibits 6 as the potential energy referred to state II under the conditions of this 
motion —no change of body and surface force except at fixed supports. 

Stability of state II implies an immediate decrease in the kinetic energy following 
the projection from state II, and therefore an immediate increase of 6’. Thus stability 
means that 6V as given by (24) is positive for arbitrary éu, . If it is given that state II 
is not unstable, (24) is not negative. 

We may take the u, to be zero, as a special case, state II then being the same as 
state 1. Evidently the equilibrium in state I, under the initial stress S,; , will be unstable 
when the right hand side of (24) is negative for any éu, which vanish at fixed supports. 

Che value which V, as given by (21), takes when the u, are the correct displacements 
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of state II can be reduced to a simpler and useful form. Writing —/ for the last integral 
in (21) we have 


a 


] “~y 
—J] = 5 | S5.U; Ui k dw 


am 


8 or a i 

=F | (S;,U;..U:),; do — 4 CS; ,t;,%),;U; dw 
a aa i ; 

=> S;.U;,uw; da — 5 (S;,;.4),;u; dw. 


Using (11) and (10) respectively in the first and second of these integrals we find 


Be ie i? i 
| 7, ;uv; do + af AF uu; dw + =f 7i;,;U; dw (29) 


= 


: l 
| AT zu; do — 9 


vz Zz 


—f = 


Nie 


If in the last integral we write 


Ti7,5Us = (TU) 5 Tea; 5 


the first of the two resulting integrals will, by the divergence theorem, cancel the second 
integral on the right of (29). Then, observing that, since 7;; = 7;; and e;; is given by (12), 


7 


] [ T;jU;,; dw = : | 7,;€;; dw = U(d, 


m& 


we can rewrite (29) as 


: £ : l f ; ze 
—-J[= 5 | AT; u; do + =f AF; u; dw — Ue). 


This form is valid only when wu; are the correct displacements of state II, because (10) 
and (11) have been incorporated. Returning to (21) we have the corresponding value 
of V as 
; ar i ; 
v=—-2[ atiudo—5[ AF, u, do (30) 


s 


e 


5. Elastic buckling. We have so far been concerned with two neighboring equilibrium 
states, state I and state II, the passage from state I to state II being effected by additional 
surface and body forces AJ; , AF; . When the passage is a buckling deformation, there 
will be no change in body force (e.g. gravity), but for very exceptional problems, and 
we may take AF; = 0. The surface forces may be taken to change at supports (e.g. the 
transverse reactions induced when buckling occurs in a column with one end clamped, 
the other pinned), but not elsewhere (the loads remain unchanged during buckling, 
moving with the particles they act on). 

Let the supports be such that no work is done by the reactions on the buckling dis- 
placements (as is true for the common boundary conditions of bars and plates. If work 
is done, as by elastic restraining moments, the elastic restraints may be included in the 
structure, and thezr fixed supports are then of the assumed type). 


GENERAL THEORY OF ELASTIC STABILITY 379 


195 


Then each integral in (30) vanishes, and therefore V = 0. Thus (21) now yields 
V = U0 +5] Sit.mir do = 0 (31) 


which is a generalized energy relation valid when S;, is a critical state of stress, and the 
u, are actual buckling displacements. 

But U’(e) is necessarily positive for any e,; since it has the form of the strain energy 
in the absence of initial stress. We have therefore from (31) 


“en | Sis, Ui. dw < 0 (32) 


and hence for a critical state I is positive. 
The initial stress S;, , now of course a critical state of stress under which buckling 


from state I to state II is possible, can be represented as 


SS, = TS; (33) 


ik 9 
where S;, is a non-critical stress tensor having the same distribution but a non-critical 
magnitude, and T is a positive multiplier. The S}, being chosen, we inquire what value 
of T corresponds to a critical state. It now follows from (32) that 


= -5 [ Si.u;, Us. dw > 0 (34) 


which defines J”. 
Introducing (33) in (31), and using the equality in (34) we have 


P= Se (35) 
I 
This holds, giving the critical value of T, when U(e) and J” are evaluated from the 


correct buckling displacements wu; . 
We now write I’ for the quantity which is calculated from the formula (35) using 
functions wu! other than wu, in U(e) and J°, leaving Sj, in the latter unchanged. We can 





then inquire what choice of u/ will yield the least value of I’. Let uf = u; + du; , u, 
being the correct buckling displacements, and correspondingly write [’ = T + 6r. Then 
‘le U_ U I° + (6U — Tel’ 
p+ er = UO+ sl _ VO + Te 4 rn. (36) 
I+ $l r+ él 


Let the variations du; satisfy the same boundary conditions as the u; . Then 


| AT’; 6u; do = 0, 
and we have from (22) 


6V = 6U + Six 6(u;, Us.) deo 


Nie 


a 0 
aU + 50 [Sh ilus.jus.s) de 


= 65U -T 6. 
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by the definition of J” in (34). Putting this in the last member of (36), and at the same 
time replacing U’(e) by TJ” from (35), we find 
a 6 V ~~ 
r+ or 
In section 4 we found that when state II is itself stable 6V' as given by (24) is positive. 
We now define state II as a buckled state which is itself stable with respect to dis- 
turbances which project it into a different configuration (the u; + 6u, not merely pro- 
portional to the u;). These disturbances correspond to the ‘‘non-sinusoidal disturbances’’ 
of the column in section 1. Our “incorrect”’ displacements u; + 6u; are of this character, 
and therefore 6V is positive. Since, by (24), 6V is zero in the first order (in éu,) and 
positive in the second, we have from (37) that 6I is also zero in the first order, positive 
in the second order, and hence that the correct value [ is the lower bound of the ap- 
proximations I’ 
Unlike 6V, 6© does not terminate with the terms of the second order. It is con- 
so that the 


ceivable that the é6u; could be chosen (not small compared with the u 
denominator in (37) becomes negative. Then the approximation I’ would be smaller 
than I. 
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Summary. Earlier results {1,2]° on safe loads for a Prandtl-Reuss material subject to 
surface tractions or displacements which increase in ratio are here extended to any 
perfectly plastic material and any history of loading. 

1. Introduction. The general limit design problem is concerned with a body or 
assemblage of bodies made of perfectly plastic (i.e. non-workhardening) material and 
subject to an arbitrary history of loading. In many cases, only the extreme values of 
the loads are given, but the order in which these loads are applied to the body is not 
specified. An important question is whether the body will ‘‘collapse’’, that is deform 
ippreciably under essentially constant loads, or whether its deformation will be con- 
tained although substantial portions may go plastic. In engineering design, the actual 
problem is to insure a reasonable margin of safety against such collapse. 

In the present paper a somewhat restricted form of this general problem is discussed: 
the actual history of loading is assumed to be completely specified rather than only 
the extreme values of all loads. The given loading history may be very simple; for 
instance, all loads may increase so as to preserve their ratios (proportional loading). 
On the other hand, the loading program may be very elaborate; additional loads may 
be superimposed on a state of initial stress, as is the case when traffic loads come on a 


bridge, or when torsion and bending are applied to a shot peened and hence prestressed 
axle 

The boundary conditions are assumed to be of the stress type for a single body or 
issemblage of bodies. At each point of the surface of the body or assemblage of bodies 
each component of the surface traction is specified except when the corresponding 
component of the displacement is prescribed to be zero. 

Several terms must be introduced, and a number of concepts must be discussed 
before the main theorems can be stated and proved. 

2. Perfect plasticity. In the following discussion, the general stress-strain relation for 
a perfectly plastic material will be used so that the results will apply to a wide variety 
of materials. By definition, a perfectly plastic material in simple tension has a stress- 
strain diagram of the form shown in Fig. 1. The essential feature of this diagram is the 
flat yield which produces a sharp boundary between elastic behavior and unrestricted 
plastic flow at points such as B in Fig. 1. 

To describe this kind of behavior mathematically for more general types of loading, 
it is convenient to use tensor notation. Latin subscripts take the range of values 1, 2, 3, 
‘Received April 4, 1951. The results presented in this paper were obtained in the course of research 
conducted under Contract N7onr-358 between the Office of Naval Research and Brown University. 

Numbers in square brackets refer to the Bibliography at the end of the paper. 
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and the summation convention concerning repeated letter subscripts is adopted. The 
coordinates x; used in the following are rectangular and Cartesian. Differentiation with 
respect to a coordinate is indicated by a comma followed by the appropriate subscript 
(u;,; = Ou;/O02;). 

The mechanical behavior of a perfectly plastic material is completely characterized 
by its yield function. For a homogeneous material, the yield function f depends only 
on the nine stress components o;; ; it is positive definite and is symmetric with respect 
to the conjugate shearing stresses o;; and o;; (¢ * j) which are formally treated as inde- 


pendent variables. Plastic flow can occur only under states of stress for which f = 1. 


o} 
A B 


f 


———____+ 





FIG. | 


rE 





States of stress for which f > 1 are not possible in a perfectly plastic material. For the 


completely stress-free state f = 0; for any other state of stress within the elastic range 
0 < f < 1. In the following, states of stress for which f < 1 will be called safe. 

The most frequently used form of the yield function is f = s,;s;;/2k’, where s,,; = 
o;; — o,,6,; is the stress deviation and k the yield stress in simple shear. However, more 


complicated forms may be used to represent, for instance, various types of anisotropy. 
For a non-homogeneous material, the yield function may vary from particle to particle. 

The stress-strain law of a perfectly plastic material does not contain the strain itself 
but only the strain rate. Since viscosity effects are disregarded, this law must contain 
the rates of stress and strain in a homogeneous manner. The strain rate e,, can be 


decomposed into an elastic component ¢{; and a plastic component ¢€7; : 
© Dp 
Cee Sag + €ij - (1) 


The elastic strain rate ef; is related to the rate of stress o{; by the generalized form 
of Hooke’s law. For the following, it is not necessary to write down this relation between 
e;; and o/; ; it suffices to note that 

€;;0;; > 0 except when i; = oF; = 0. (2) 

The plastic strain rate is related to the state of stress. Since viscosity effects are neg- 
lected, the stress components must be homogeneous of the order zero in the components 
of the plastic strain rate. In other terms, the stress tensor o,;; determines the tensor of 
the plastic strain rate «2; only to within an arbitrary factor. 

The following geometric representation of states of stress and plastic strain rate is 
often useful. In a nine-dimensional Euclidean space, consider a fixed system of rect- 
angular Cartesian coordinates. The state of stress o,;; will be represented by the point 
with coordinates proportional to o;; . The plastic strain rate ¢{; will be represented by 
the ray with direction cosines proportional to ¢;; . This manner of representing the 
plastic strain rate by a ray rather than a point is suggested by the fact that the state 
of stress determines the plastic strain rate only to within an arbitrary factor. 
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The yield condition f = 1 defines a surface in this nine-dimensional space. We assume 
this yield surface to be convex and to possess a unique normal in each of its points. 
The mechanical significance of this assumption will become clear from the following 
discussion. 

Consider a state of stress for which f = 1. As has been shown in earlier papers [3, 4], 
any plastic strain rate ¢€; associated with this state of stress is represented by the ray 
which has the direction of the exterior normal of the yield surface at the point o;; . 
From this fact and the assumption regarding the yield surface, there follow two im- 
portant lemmas. 


Lemma 1. The stress rate o{; and the plastic strain rate ¢, satisfy 
oie; = 0. (3) 


Proof. If the plastic strain rate is not to vanish, the stress rate must correspond to 
the change from one point on the yield surface to a neighboring point. Thus, the stress 
rate is represented by a vector which is tangential to the yield surface, whereas the plastic 
strain rate is represented by a ray normal to the yield surface. Equation (3) expresses 
the orthogonality between this vector and ray. The special case o/; = 0 also satisfies 
Kg. (3). 

Lemma 2. There exists a function F, homogeneous of the first order in the com- 
ponents of the plastic strain rate, such that in flow with the plastic strain rate €; , 
energy is dissipated at the rate F'(e;;), whereas, for any safe state of stress oj; , 


‘7 


oii€:3 < F(é;). (4) 

Proof. lf the plastic strain rate ¢?; is associated with the stress o;; , the rate of dissipa- 

tion of energy is o,;€; . Since the stress components are homogeneous of the order zero 
II 


FIG. 2 





in the components of the plastic strain rate, the dissipated energy is homogeneous of 
the first order in these strain rate components. Let IT be the tangent plane of the yield 
surface at the point o,; and h its distance from the origin 0 (Fig. 2). The rate of dis- 
sipation of energy is then given by (e?,€7;)'”"h. The inequality (4) expresses the fact 
that any safe state of stress is represented by a point which lies on the same side of II 
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as the origin. (Note that the function F is the supporting function of the convex yield 


surface 

3. Collapse. From the point of view of the practical engineer, the term collapse 
implies that appreciable changes in the geometry of a structure will oceur under es- 
senti:lly constant loads. For the purpose of the following discussion, however, it is 
more convenient to use the term collapse to refer to conditions for which plastic flow 
would occur under constant loads if the accompanying change in the geometry of the struc- 


} 


ture or body were egarded. In the discussion of this type of collapse, the equilibrium 


conditions can be set up for the undeformed rather than the deformed body. This 
obviously represents a great simplification because the deformation occurring during 
collapse is not known beforehand but constitutes one of the unknowns of the collapse 
problem 

To illustrate the relation between these two definitions of collapse, consider first a 
hollow sphere of uniform wall thickness under gradually increasing interior pressure 
The material at the interior surface reaches the yield limit first, but its plastic deforma- 
tion is contained by the surrounding shell of still elastic material. As the pressure con- 
tinues to increase, the boundary between the elastic and plastic regions moves tawards 
the exterior surface; only when it has reached this surface does large plastic flow become 
possible. Let us now compare the simplified stress analysis which neglects all changes of 
geometry during the loading process to the complete analysis which takes account of 
these changes. Even in the elastic range, and also in the subsequent range ol contained 
plastic deformation, the hollow sphere expands and its wall thickness diminishes some- 
what. If those effects are taken into account, the pressure for which the entire sphere 
becomes plastic is found somewhat smaller than when they are neglected. However, the 
difference between the two pressure values is small because the deformations occurring 
up to the end of the range of contained plastic deformation are of the order of magnitude 
of elastic deformations. Continuing the analysis into the subsequent range of unrestricted 
plastic flow, we find that the simplified analysis predicts flow under constant interior 
complete analysis predicts flow under gradually decreasing pres- 


pressure, whereas the 
as this 


sure. In either case the sphere loses its usefulness as a pressure vessel. As fat 
first example is concerned, the two definitions of collapse are therefore in substantial 
agreement as to the pressure under which collapse sets in. 

Next, consider a blunt, rigid wedge which is pressed against the flat surface of a 
large block of perfectly plastic material. It is obvious that the term “collapse” can not 
be applied to this problem in the sense attributed to this term by the practical engineer. 
Indeed, us the wedge is pressed into the plastic material, the area of contact continues 
to increase; thus, the load on the wedge must be increased steadily if plastic flow is to 
be maintained. On the other hand, there exists a well-defined “collapse load’’ under 
which plastic flow would continue if all changes in geometry could be disregarded. The 
physical meaning of this ‘collapse load”’ is less clear than that of the collapse pressure 
obtained in the preceding example. Obviously, the collapse load in the case of the wedge 
is not the load under which the first indentation occurs, because some identation takes 
place even in the elastic range. It is likely, however, that the collapse load indicates 
the load intensity at which the permanent indentation begins to increase considerably 
faster than the elastic identation. Thus, the collapse load in the case of the wedge has 
the nature of a conventional yield limit, whereas the collapse pressure in the case of the 


hollow sphere represents a natural yield limit. 
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To arrive at a mathematical characterization of collapse, let the velocities be de- 
noted by v, , the surface tractions by 7’; , and the body forces (per unit volume) by 
fF, . Furthermore, use the prime to indicate rates of change and the superscript c to 
refer to collapse. From the definition of collapse introduced above it follows then that 


during collapse 


| T’v, dS = 0 and Fr = @ for some v; ~ 0, (5) 


where the integration is extended over the surface S which bounds the considered body 
or assemblage of bodies. Indeed, since collapse is to occur under constant loads, the 
rates F%° must vanish throughout V; moreover, for those components of the surface 
traction which are prescribed by the boundary conditions the rates 7% must vanish, 
whereas the surface velocities v{ corresponding to the remaining components of the 
surface traction must vanish according to our definition of stress boundary conditions. 
Thus, the integral in the first Eq. (5) is seen to vanish. 

An expression for the rate at which work is done during collapse is useful. If the 
velocities v{ considered as functions of the coordinates are continuous and have con- 


tinuous first derivatives, the principle of virtual displacements yields the following 


equation: 
| Tvias+ | Pvav = | ote, aV. (6) 


4. Admissible states. Consider first a state of stress for which the components 
are continuous functions of the coordinates. Such a state is called statically admissible 


og 


if it satisfies (7) the conditions of equilibrium 


Ci tf; = 0 (7) 
throughout V and (77) the boundary condition 
a,;n; = T,; S) 


on those portions of the surface where the component 7’, of the surface traction is given. 
In (8), the unit vector along the exterior normal of S is denoted by n, . 

The preceding definition may be generalized to include stress fields with a finite 
number of surfaces of discontinuity. On either side of such a surface, the stresses must 
then satisfy (7). Moreover, if n* denotes the unit normal vector of the surface of dis- 
continuity, the expression ¢,,n* must have the same value whether it is evaluated from 
the stresses on one or the other side of the surface of discontinuity. 

A velocity field v, is called kinematically admissible if the velocity component v, 
vanishes on those portions of the surface S where the corresponding component 7’; of 
the surface traction is not prescribed. A kinematically admissible velocity field may 
represent rigid body motions for certain portions of the body and genuine deformations 
for the remainder. Special discontinuous velocity fields are also permissible and often 
useful; they will be discussed in some detail later. For the present, however, we consider 
only kinematically admissible velocity fields for which the velocity components are 
continuous functions of the coordinates. 

Such a velocity field is said to define a kinematically admissible state of collapse if the 
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rate at which the actual surface tractions and body forces do work on the velocities v; 
equals or exceeds the rate of dissipation of energy computed from the strain rates 


esi = 405.5 + 05,,) (9) 


treated as purely plastic strain rates. Thus, for a kinematically admissible collapse state 


a . 


[ rotas+ | Fotav > | Feéyav. (10) 


© d 


5. Collapse theorems. The following theorems were previously established for special 
cases [1, 2]; they can now be shown to hold generally. 


Theorem 1. If all changes in geometry occurring during collapse are neglected, all 
stresses are found to remain constant during collapse. 


Proof. Applying the principle of virtual work to the velocity field and the rates of 
change of the surface tractions, body forces, stresses and strains during collpase, we 
obtain 


a . . 


| Ti dS + | FrvidV = | ofes, av. (11) 
According to (5), the left-hand side of (11) must vanish for the considered collapse 
state. The strain rate on the right-hand side of (11) can be decomposed into its elastic 
and plastic components; thus, 

| otiess dV + | offs aV = 0. (12) 
The second integral in (12) vanishes according to (3). The relation (2) shows then that 
vanishes throughout V. 


ye 


(12) can be satisfied only if the stress rate of; 


Theorem 2. If a safe statically admissible state of stress can be found at each stage 


of loading, collapse will not occur under the given loading schedule. 


Proof. Suppose this theorem to be false. Then, at some stage of loading, a collapse 
state v{ would exist although a safe statically admissible state of stress o;; could be found. 
Applying the principle of virtual work to the actual surface tractions 7; and body forces 
F* at this collapse stage, the stresses o;; with which these are in equilibrium, the velocities 
v; , and the corresponding strain rates ¢«;; , we obtain 


+ 


Tv dS + / Fv dV = | o°,¢¢, aV. (13) 


ny ny 


According to Theorem 1, the stresses and hence the elastic strain remain constant 
during the collapse described by the field v; . Thus, the strain rate ¢{; is purely plastic. 
From the first part of Lemma 2 it follows then that the right-hand side of (6) can be 
written as f F(e;;) dV. Combining this form of (6) with (13), we obtain 


| oi, dV = | F(é,) dV (14) 
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which is in contradiction to the second part of Lemma 2, since the strain rate ¢{; is 
purely plastic. 


Theorem 3. As long as collapse does not occur, a safe statically admissible state of 
stress can be found at each stage of loading. 

To prove this converse of Theorem 2, consider a generic stage of loading defined by 
the surface tractions 7’; and the body forces F; . If the actual stresses at this stage of 
loading are denoted by a,; , we have f(o;;) < 1, because collapse is not supposed to occur. 

Consider now the loading specified by V7’; and NF; , where N is constant throughout 
the body. For NV = 1, collapse does not occur according to the condition of the theorem. 
Since the equations of equilibrium are linear in the stresses, body forces, and surface 
tractions, the stresses No,; will be in equilibrium with N7; and NF; . Moreover, it 
follows from the convexity of the yield surface that f(No;;) < lif N < 1. Thus, collapse 
can occur under the loads NT’; , NF; only if N > 1. Let o{; denote the actual stresses 
for such a state of collapse. These stresses are in equilibrium with N7’; , NF; and satisfy 
f(o°;) < 1. Therefore, the stresses o{;/N are in equilibrium with 7’; , F; and satisfy 
f(o°,/N) < 1; in other terms, these stresses define a safe statically admissible state of 
stress for the loads 7’; , F; . 


Theorem 4. If a kinematically admissible collapse state can be found at any stage of 
loading, collapse must impend or have taken place previously. 


Proof. Suppose this theorem to be false. According to Theorem 3, a safe state of 
stress o:, could then be found in spite of the existence of a kinematically admissible 
collapse state v; . Applying the principle of virtual work to the actual surface tractions 
7, and body forces Ff, at the considered stage of loading, the stresses o;; with which 
these are in equilibrium, the velocities v; and the corresponding strain rates e;; , we 


obtain 
| Twi dS+ | Fw: dV = [oie dV. (15) 


On the other hand, we may use (10) since v; is a kinematically admissible collapse state. 
Thus, 

| a',¢,,;dV > [ F(é;) dV (16) 
which is in contradiction to (4) since the strain rates ¢;; of a kinematically admissible 
collapse state are treated as purely plastic. 

6. Discontinuous velocity fields. It is often useful to consider discontinuous velocity 
fields. However, it should be kept in mind that in plastic flow, as distinct from fracture, 
actual discontinuities cannot occur across a fixed surface. The type of discontinuity to 
be considered in the following is simply an idealization of a continuous distribution in 
which the velocity changes very rapidly across a thin transition layer (Fig. 3). This 
idealization is permissible provided the stresses on the assumed discontinuity surface 
are chosen as the limiting values of the stresses on the surfaces bounding the transition 
layer as the thickness of this layer approaches zero. For plane plastic flow in a Prandtl- 

teuss material, for instance, the line of discontinuity must be a shear line for each of 
the stress fields on the two sides of the line of discontinuity. If the yield function f 
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depends upon the mean normal stress, it will be found that a discontinuity in tangential 


velocity implies separation or overlap of the material on the two sides of the discontinuity. 


In such a case, the actual transition layer must have appreciable thickness, but the 


scontinuity surface may still be useful.* 


idealization to a d 
The theorems of Sec. 5 are obviously valid in the presence of a transition layer. They 


will, therefore, remain valid in the limit as the thickness of the transition layer approaches 


FIG. 3 





zero, provided it is kept in mind that the rate of dissipation of energy in the transition 
layer approaches a finite value in the limit. If the transition layer is replaced by a surface 
of discontinuity, the expression | F(e;;) dV must, therefore, be replaced by 


| F(e,,) dV + | Tdv, dA (17) 
where dA is the element of area of the discontinuity surface, 7; the traction and Av, 
the velocity jump across this surface. Thus, in the definition (10) of a kinematically 
admissible collapse state, the right-hand side must be replaced by (17). 

Actual discontinuity of velocity can occur in the case of an assemblage of bodies; 
for instance slip may occur between a punch and the indented material. If there is no 
friction between the bodies of an assemblage, Theorems | through 4 remain valid in spite 
of such discontinuities in the velocity because no energy is dissipated on the contact 
surface in the absence of friction. 

7. Additional theorems. The following theorems are intuitively obvious but their 
statement and indication of proof seems worthwhile. 


Theorem 5. Addition to the body of (weightless) material cannot result in a lower 
collapse load. The proof follows directly from the fact that the collapse stresses a; 
for the original body and zero stresses in the added material constitute a statically 


admissible state for the new body. 


Corollary. Removal of material cannot increase the collapse load. 
If the yield surface of one material contains that of a second material, the first 
material will be said to have higher yield strength than the second. 


Theorem 6. Increasing the yield strength of the material in any region of a perfectly 
plastic body cannot weaken the body. 

The proof follows from the fact that any statically admissible state of stress which 
is safe for the unimproved body is also safe for the improved body. 

*Application to soil mechanics provides an excellent illustration and will be treated in a separate 


paper 
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Theorem 7. Initial stresses or deformations have no effect on collapse providing the 
geometry is essentially unaltered. 

To prove this, we note that an initial or residual state does not affect any of the 
equations or statements made. This means, for example, that settlement of supports of 
L continuous structure or initial plastic torsion of a bar subsequently bent or pulled 
does not affect the limit loads provided the geometry is not changed appreciably. 

It is probably best to end on a note of caution. Just as in elasticity, the concept of 


an essentially unaltered geometry rules out buckling which must, therefore, be studied 


separately. 
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AN INTEGRAL EQUATION APPROACH TO THE PROBLEM OF 
WAVE PROPAGATION OVER AN IRREGULAR SURFACE* 


BY 
GEORGE A, HUFFORD 
National Bureau of Standards 


In this paper we propose to outline a certain approximate technique for solving the 
wave equation under the following conditions. Consider a surface S which resembles 
somewhat the surface drawn in Fig. 1; it stretches out to infinity along a horizontal 








Fic. 1. 


plane and is composed of “‘large’’ corrugations or bumps—large in the sense that the 
radius of curvature at any point of S is much larger than a wavelength. Above S there 
is a field y (x, y, ze ‘** which varies harmonically in time with the radial frequency w 
and which satisfies the wave equation 


Viv + ky = —4ar, (1) 


k being the propagation constant w/c and + the distribution of sources. Finally, at the 
surface S the field satisfies a homogeneous boundary condition of the form 


oy _ 
oe = —ik by, (2) 


where 6 is a certain (complex) constant of proportionality. 


*Received Nov. 22, 1950. 
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As we have formulated it, this problem first arose as a simplification of the problem 
of radio wave propagation over that irregular and inhomogeneous interface that is our 
Earth. At the higher frequencies being used today, radio wave transmission is greatly 
affected by these irregularities, and it is becoming increasingly important to find some 
way of estimating the resulting field strengths. 

The simplifications in the present case involve (1) a scalar wave phenomenon, and 
2) a homogeneous boundary condition. But, while the remaining problem is not itself 
without merit, it might be remarked that both assumptions can be looked upon as 
reasonable approximations to the original radio wave problem. The first follows trom 
the fact that a vertically polarized source will give rise in the main to a vertically 
polarized field, and a horizontally polarized source to a horizontally polarized field. 
Thus the vertical electric field or the vertical magnetic field or a Hertz potential, all of 
which satisfy the scalar wave equation, will serve to characterize the field. As for the 
second assumption, divers arguments that lend credence to it have been presented by 
Schelkunoff [1], Leontovich [2], Leontovich and Fock [3], Feinberg [4], and the present 


‘ 


author [5]. These show that in a vertically polarized field something equivalent to Eq. (2 
will be satisfied with 

. = oe | ‘a 

= , 2) 


n n’ 
where 7 is the complex dielectric constant of the earth 


€ e Cx 
j==—4+ t=, 
€; €yW 


being the dielectric constant and o the conductivity) and vn is the coefficient of 


€/€ 
refraction equal to n' *. Similarly in a horizontally polarized field Eq. (2) will be satisfied 
with 


6= (yn — 1)" =n. 


Previous authors in their attacks on the wave equation in the presence of irregular 
boundaries have for the most part used a perturbation method. This method was 
originated by Brillouin [6] and later developed by Feshbach [7] and Cabrera {8}. It 
begins by deriving an integral equation and then solves this by approximating the char- 
acteristic functions through known solutions to problems involving regular boundaries. 
However, it was pointed out by Feinberg [4] that in the special case of radio transmission 
this same integral equation could be reduced to a much simpler form. Subsequently, 
the present author in his Master’s thesis [5] tried to simplify the integral equation 
further in the hopes that it could then be directly solved by numerical methods. The 
present paper is a recapitulation and a continuation of this approach. 

The integral equation. According to Green’s theorem, if ¢ and y are any two functions 
continuous throughout a volume V, then [9] 

[ure-o'nw=[(o-vu 
where S* is the surface bounding V. Now let us consider the volume which is bounded 
as in Fig. 1 by the surface S and the large hemisphere-like surface >. Let P be some 
fixed point on S and Q the variable point of integration of either the volume integral or 
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the surface integral. Then suppose that ¥(Q) satisfies Eqs. (1) and (2) while ¢(Q) is 
the Green’s function e'""*/r, , r. being the radial distance PQ. Of course, since ¢(Q) as 
thus defined has a singularity at P, we must exclude P from the volume V by indenting 
it with a very small hemisphere co. 

Under these hypotheses the volume integral of Eq. (5) becomes 


ikr 


ir | 7(Q) —— dv = 4ny,(P), (6) 
> 


an 2 


say. This function y,(/) is the field that would be obtained at P if there were no earth 
to contend with. It is the ‘free-space field”? due to the source distribution r. 

The surface S* consists of three parts: 2, ¢, and S. Since y must satisfy the radiation 
condition [10], it is clear that as the radius of = becomes larger and larger, the integral 
over = vanishes. If the radius of o tends to zero, then because of the singularity in (Q) 
the integral over ¢ approaches 27y(P). Finally we know that Eq. (2) is satisfied at all 
points on S. Thus, after these considerations and after a little rearrangement, Eq. (5) 


hecomes 11] 
F s ii, 74 
) — on (P) 2 ia ( J ) a | ” 
VP) = 2(P) + 5 I ¥Q) — E +(1+5-)5,| (7) 


This as you see, is an integral equation which defines uniquely the unknown y(P) at 
all points on S. We want now to simplify it with a few judicious approximations. 
First, we suppose that the sources are all located on some antenna structure which 
is erected at a point O on S. Then we write 
tkro 
vo(P) = g(P)—, (8) 
0 
where r, is the radial distance OP. The function g(P) is, in a way, the antenna pattern; 
it represents the gain of the antenna over an isotropic radiator at 0. 
We also introduce an attentuation function W(P) which is defined so that 
ikr 


v(P) = 2p) —. (9) 


Substituting these two equations into Eq. (7) we have 


wir) = oP) + f wager [54 (14 2) | a (10) 
2r Js TNs kr./ an 

Here now is the crux of the matter. Because of the factor exp ik(7, + r. — ro), the 
integrand in this last equation oscillates very rapidly from one point on S to another. 
This fact implies that we might use Kelvin’s principle of stationary phases; that is, we 
could approximate the integral by summing up only those contributions which come 
neighborhoods of points where the phase of the integrand is stationary. These 
points in general form a discrete set which are exactly the “points of reflection’’ used in 
geometrical optics, and indeed, doing this should give us precisely the geometrical optics 
solution. It may be remarked that this very solution has already been used with some 
success by MecPetrie and Ford [12] and by Shelleng, Burrows, and Ferrell [13] to analyze 
some actual data, while Keller and Keller [14] have devised several formulas to be used 


from the 


in this connection. 
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However, we cannot in our problem use this particular approximation, for when the 
geometrical rays are at nearly grazing angles of incidence or, even worse, when we must 
consider diffraction regions, then it is not nearly accurate enough. The reason for this 
is that at points between O and P the phase, while not necessarily stationary, is still 
but slowly varying. Thus to improve upon the geometrical optics solution we must 
reduce the integral in Eq. (10) not to the sum of discrete contributions, but to a line 
integral from O to P. 

To do this we first project everything onto a horizontal plane S’. Then the integral 


can be written in the form 


. Peo ers 
T= | FQ) tS, 
Js: ri 


where the primes are used to represent horizontal projections. Now let us construct on 
S’ a set of elliptic coordinates (u, v) defined by 

r,coshu=r34+7t, 7, COSY = 7, — 17, 
The differential area is r3r{ du dv and the integral becomes 


a «x a* 


[= | exp [zkrj(cosh u — 1)] du | F(u, v) dv. 


J “0 


This, as an integral in u, has the phase kr§ (cosh wu — 1) which is stationary only at 
the point u = 0. Immediately, therefore, we can write down the approximation [15] 


ee ee 
kro Jo 

This approximation is good, of course, only if [> F(u, v) dv is a slowly varying function 
of wu. But this will always be the case if S is smooth enough. If, on the other hand, there 
are projections on S which are away from the line OP but which nevertheless contain 
points of reflection, then surely these must also be taken into account in evaluating the 
integral. In what follows we shall assume that no such projections exist. 

The line u = 0,0 < v < az, to which we have reduced the surface integral, is the line 
segment O’P’. Adopting at this point a more suitable one dimensional notation, let us 
denote the distance rj by the letter x and the distance of the point (0, v) from O by s. 
Then, 

s = 32(1 — cos»), 


so that our approximation becomes 
Qn \'/? aah: as ds 
eel ee 
ka Jo [s(x — s)] 

The final form of our integral equation can be further simplified by making several 
more approximations in F(s). Thus we assume that da/da’ is close to unity, that 1/kr. 
dr,/dn is negligible, and that whenever 7, , 7; , 72 appear in the modulus of F(s) they can 
be replaced by the horizontal distances x, s, x — s, respectively. Equation (10) then 
becomes 


y am \ — ,-in/4 zy" Po E »( #1) 2 roreensl 1M . 
W(x) = g(x) — « (# | W(s)\ 6 + an) pean ds. (11) 


0 
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This is the integral equation we are proposing. It is, so to speak, an approximate 
representation of our problem which has the advantage of simplicity while still retaining 
most of the characteristics upon which one intuitively feels the problem should depend. 
Actually, the reduction to one dimension is not an unfamiliar concept. Experimental 
and theoretical investigators have for a long time drawn profiles of the earth and at- 
tempted to correlate the profile with measured field strengths. Following our own 





























Fia. 2. 


proposal, we would draw a profile similar to the one in Fig. 2; compute from it the 
quantities 7) , 7: , T2 , Or2/8n; substitute these in Eq. (11); and then, by using suitable 
numerical methods, solve for W(z). 

Indeed, let us represent the profile in question by the function ¢(x). This function 
will be the elevation of the earth above a datum plane passing through the point O. 
Then we can write 

2 2 1/2 (x) 
mn=(r+°@)" =zr+ “Oe 
and similar expressions for r,; and r, . A little algebraic manipulation of these expressions 


gives us 


eer me! pe. iy 
r+ rs r= tee | of (12) 
Similarly, 

nt th nt ees a 

an = Sin (8 — a) = (8) ae (13) 
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The field above the earth. We have thus found a procedure for evaluating W(x) and 
hence the field y at all points of the earth’s surface. It remains for us to show how this 
can be used to find the field at points in the space above this surface. 

Using Green’s theorem again together with Eq. (1) we find 


; ; pt gh v1) 
T T\ di — 7 | | si 7 — . 
' | @ R aha I. | HQ) on R R, oan 


where this time P is in the space above S. The quantity FR, is the radial distance PQ; we 
use the capital letter to stress the fact that P has been raised off the earth. Note that to 
exclude P from the integral of Green’s theorem we now must construct a complete sphere 
about P rather than a hemisphere. It is for this reason that the 2m(P) has changed to 
tory (P). 

Now this equation, too, contains a surface integral whose integrand oscillates rapidly 
everywhere except possibly between O and P. Hence we can follow exactly the same 
steps we used before to change the integral to an approximate one-dimensional integral. 
If we replace the volume integral by ¥(P) and use Eq. (2) to dispose of dy/dn, we shall 


1/2 
? x 
Rs | | ds. (14) 
eae 


and dR,/dn can be approximated by suitable modifications of 


finally arrive at the equation 


ae eo y a” y(s a) 
er oe (¢ | rine + oP 


“0 \ 


The quantities 7, + R 
Eqs. (12) and (13). 

A plane earth. Although Eqs. (11) and (14) were derived largely to provide numerical 
solutions, it is perhaps of some interest to examine a few special cases which can be solved 
analytically. 

One such case is that of a plane, homogeneous earth with an isotropic antenna at a 


height h, above the point O. Then dr,/dn = 0,7, + r2 — 7% = O, and 


Vola) = (0* + Ri)? exp [ak(a® + hi)'] = = exp n(x + x)| 
x or 


or 

g(x) = exp (ikh;/2z). (15) 
Now it may be argued that this last approximation can hardly be valid since at the 
point O where x = 0, the error becomes infinite. But while this will certainly affect the 
solution near O it will actually have very little affect further away. This is because the 
phase of the error also becomes infinite causing the integrated error to be small. 


Introducing, then, these quantities into Eq. (11) we have 


. _ Pye k 1/2 az : z 1/2 
W(a) = exp (tkh,/2x) — e*’ (J ) | WW (| mes | ds. (16) 
2 Jo as 2 


This we can simplify further with the use of Sommerfeld’s complex numerical distance. 


We define 
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whereupon Eq. (16) becomes 


—" 1/2 
Wp) =e °°? + ig” | wo] 5 p | dy, (17) 
0 vip — v 


vhere a kh,6é 

by using the 
Liouville-Neumann series or by a method similar to that used in solving Abel’s integral 
equation. We shall find it most convenient to use Laplace transforms. This we do first 


This integral equation can now be solved in any of several ways 


by defining 


U(p) = p '*W(p); 1S) 
then Kq. (17) becomes 
— ; i U(v) 
U(p) = p —" + ir ash | 5 Gv. (19) 
F ja OOD 


In this equation the integral is of the Faltung or convolution type, so that when we apply 
the Laplace transformation to both sides we obtain [16] 
a\'’* 
L{U} = ( ew” +. LIU} 
Pp 


or, solving for L{U} 


L{U) =< =, = (*) ce" 4 (*) ———. (20) 
mt a p p - =% 


It only remains for us to find the inverse transform of this function. But the first 
term is the transform of p '*e~*’’, while the second term we can write as i(#/p)'*/(p' *) 
which is the Laplace transform of [17] 


ip”? | e' V(t) dt, 


“oO 


where 


joy, Dee i <a, 
' 
\p- i ) 


This last, of course, assumes that @ is real. But the functions that we obtain in the 
end will all be analytic in a and hence valid also for complex values.) Thus 


. 1/2 a*/4 . 1/2 P = t ; _— 
U(e) = p + ip | exp | -£ + 2 «| dt 


Vid) = L {fp} = L" 
“ae i> es. 


(21) 


= fe ee wr a 
=p ec 4 tg'e? erfe (—io' - =a), 
p 
where 
Qf 


erfe (x) = -773 | ee’ du. 
T 


“2 
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Or, if we define 


w= (it) = (1+), (22) 
p 6x 
we have finally 
W(p) = p'?U(p) = exp [ikh{/2x][1 + i(axp)'e~” erfe (—iw'”’)]. (23) 


To complete our analysis, suppose that we raise the receiving antenna to a height 
h, . From Eq. (14) together with the usual approximations for R, and dR,/dn, we have 


ikRo 
WP) = = 
(24) 
(2) ‘ize [ WO) tkh, _ | “1/2 *h36°/ | 
Th * ea ee rd ee exp [—k*h26"/4(p — v)] dy, 


where FR, is the radial distance from the isotropic antenna to the point P. 
In order to evaluate this integral we shall again make use of the Faltung theorem. 
Denoting the integral by J and remembering Eq. (20), 


L{I} = L{p™ Win tA(t _ thy) exp (—Kn33'/49)} 
2p 


| as 2 
= im exp [—k(h, + h.) &p'”]| =z - as I. 
I [ 1 2) Of ] p” p’ nae 
Thus the Laplace transform of J is represented as the sum of two terms, the first of 
which is the transform of i(x/p)'’? exp [—k*(h, + h)*6°/4p], while the second is exactly 
the same transform as that in Eq. (20) except for the factor —72r'’’ and except that 
h, has been replaced by h, + h, . Therefore we have immediately 


pikRo 
> a oe 
we) = Ry 
(26) 
— : exp {ik[a + (h, + h.)?/2zx]}{1 — 2[1 + 2(p)'e~” erfce (—iw)'”’]}, 
where now 
h, +h ) se 
cod ae 2 
u o(1 - om ‘ (27) 


The problem of electromagnetic radiation from a vertical Hertzian dipole over a 
plane and homogeneous earth was first successfully attacked by Sommerfeld in 1909, 
and since then many authors have treated its various aspects. To compare our Eq. (26) 
with the classical solutions, we quote here an approximation due to Norton: [18] 


ikR ikRo’ 
e 


V(P) = . ie {1 — 2[1 + i(ap)'’e” erfc (—iw) 


“x * 
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where ¥(P) is the Hertz potential at P, 
p = tkRi/2n’, 
w = p[l + n(hi + hy)/Rel’, 
n is the coefficient of refraction that we mentioned in Eq. (3), and Rj is the radial distance 
from the dipole to the image of the point P, 
Re? =z + (hi + hz)’. 


The agreement, it will be noticed, is almost exact. 











Fig. 3. 


A spherical earth. As a last example let us consider an isotropic antenna at the surface 
of a spherical, homogeneous earth. If the sphere is of radius a, then from Fig. 3 we have 
2 


2\1/2 x (28) 


—-a>=- = 


t(z) = (a — 2°’) ~ “2, 


Substituting this approximation into Eqs. (12) and (13) gives us 


SX\TI — 
n+ —-% = ee 
2-8 
dn —- 2a = 


Now we introduce the so-called natural units &, ¢, and y, which are defined by 


t= ra**(k/2n)”, 


f= sa~*’*(k/2n)'*, (30) 


1 = i/(ka)"”, 
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and further make the substitution 
Wo) i ~* U(é). 31) 
Then our problem of a spherical earth becomes immediately one of solving the integral 


equation 


[ (x 
U® =F '" exp| (a) | 


29 
32) 
[ct np eap | (le -¢ : : (é ) HE at 
— \ ; ¢ ~ ¢) ‘ I 7 ‘ c \< _ 4 ¢ aad 4 S- 
‘ a 12 ws (29 2 ss 
Again we have a Faltung integral. If we denote by u(p) the Laplace transform ot 
U(£) and by v(p) the transform of the function & '~¢ ~*~ then taking the Laplace 
transformation of both sides of Eq. (32) and making use of the Fallung theorem and 
the rule for differentiating a transform, we find easily 
V\ Pp) ») 
Np) = -/4 i . >) 
1+. (2x) "yy u(p) — (1/2) v’(p) 
And now to find (’(£) we need only to use the complex inversion formula 
ae I [ ] 2 
(€) = a ul pye ap. 34) 
Leaving aside, for the moment, all questions of existence and convergence we remark 


that this integral can be evaluated from a knowledge of the poles of u(p). As we shall 


show later, v(p) is an integral function so that the poles of u(p) are all at the zeros of 


the denominator in Eq. (33). It follows that if these zeros are at the points p p 
n = 0,1, --: , we shall finally have 
V(p,, Je ‘ 
U(é) = > f . 35) 
, ix/47o ¥ is —— 9 am/4 ve 
 ¢ (27) yY v(p,) (1/2) v'"(p,) 


The properties of »(p). We have defined v(p) by the integral 


v(p) = | é Nak” ** dé 30) 


which, however, converges only for Re p > 0. If we make the transformation 


E =. (12/n) ek 37) 
then it follows that 
v(p) = (12/m) “ee *" "k[(12/m)'"e*" "pl, 38) 
where 
ki) = | oe -* dt. (39) 


Since this last integral is analytic for all z it defines the analytical continuation of v(p) 
to the entire p-plane. 
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Obviously A(z) is regular everywhere and is therefore an integral function. It has 
ie Taylor series expansion 

i “ae (—gi) .. I r(n/3 + 1/6) ; 

k(z , eve” Si dt == > (—z)". (40) 
— ee (in + 1) 

Now from this Taylor series, it is possible to derive an asymptotic expansion valid for 
large z. We need to know that the function I'[(z/3) + (1/6)]/T'(2 + 1) has poles at 
34 , and that everywhere to the right of the imaginary axis it has the 


mptotic expansion 


peng | TR 
rie 
/ 
4 
Wi bo 
ee, 
bo 
ww i 
a 


I 
r(2z/3 + 4/3) [ + o() 


Then it can be shown that [19] 





; 2r < < 2r 
v2 -—- argz < —, 
3 - 3 
1/2r°>_—1/2 {2/6 (93/2) ,3/2) 1/2 2r 
k(z) ~ 4 9 *[tz exp {7(2/3"")z""} + 2°] > Sargz<r, (41) 
o 
: ; ; , > 4n 
wn *li2”'”? exp {1(2/3°”)2°7} — 2" ] mr <argz< 3 
1477/3 


From this it follows that the zeros of k(z) lie on or near the radials e'**’* and e 
since it is there that the exponentials of Eq. (41) have an absolute value equal to one. 
If we return to a consideration of the original function v(p) we see from Eq. (38) that 
its zeros will lie along the radials e'**’? and e'°”’°. In fact if we let 


p9=ec "= — (42) 
then for small arguments of o and 7 we have respectively 
v(p) ~ Pea"? Lexp [i(4/32'”’)o"*] — exp (—ir/2)}, 
(43) 
v(p) ~ We exp [—1(4/39'”*)7°”?] — exp (ix/2)}, 
and further 
v'(p) ~ —2e°'" exp [i(4/32'”*)o"”’], 
v’(p) ~ —2e°'*”* exp [—7(4/39'”) 7°” ], (44) 
v''(p) ~ —4Ar Mae tN ”? exp [—1(4 30’)? “<. 


We shall need one more property of the function v(p). It can easily be shown, after 
integrating Eq. (39) by parts, that 
6k” + 22k’ +k = 0. (45) 
On the other hand 


d Ea am sare a8 «| = | p6K'” + 22k’ + k] = 0, 
dz 2 6 6 
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so that 


kk” — 5k’) “ ‘ k* = const. (46) 


The constant of integration, which can be determined by evaluating Eq. (40) and its 
derivatives at z = 0, is equal to 7/6. 
From this we can conclude that whenever k(z) = 0, then k’(z) = +2(r/3)°°. Or 


again, whenever v(p) = 0, then 
v’(p) = +27". (47) 


Furthermore, it follows from Eqs. (48) and (44) that at those zeros lying along the 


it is the positive sign that is valid, while along the radial e it is the 


32/2 


radial e’ 
negative sign. 

The field. Having thus determined the salient properties of the function v(p) we are 
now in a position to discuss the properties of u(p) and its inverse transformation. 

First of all, it can be shown from Eq. (33) and from the asymptotic expansions of 
v(p) that, except in the neighborhood of its poles, u(p) = 0(p '’*). From this it follows [20] 
that the integral of Eq. (34) converges, that it represents the function U(é), and that it 
can be expanded in the series of Eq. (35). Thus this series is established rigorously as 
the solution to Eq. (32). 

Now ordinarily y is a very small quantity so that we shall expect the zeros of the 
denominator in Eq. (33) to be very near the zeros of v(p). Thus there will be two sets 
of zeros: those in the vicinity of the radial e'**’* and those near the radial e'’*’’. And 
as a matter of fact, it follows from Eq. (47) that those zeros of v(p) lying near the radial 
* are exactly zeros of the denominator. But because they are zeros of v(p) they can 


ist 


€ 
contribute nothing to the series of Eq. (35) since the numerator of each term contains 
a factor v(p,). 

Thus we are left with only those zeros that lie near the radial « 


15a /6 


If we write 
p, = 7,e"’° and make use of the asymptotic approximations in Eqs. (43) and (44), 
then finding the zeros becomes equivalent to solving the equation 


“i(A/Bx8/2)rq8/2 +e 


x/6/ 1/6 / 1/2 
"a4 Tt, 


1 eed aie 4)! S lyri? 


nr ri” ten(=) a : 
é ae — 9 73 = —ae 1/2 (< Y 
7 (7 os) 4, YTr ” 


And now we can write down the final answer. But before we do let us note that 





(48) 


¢ = 


or alternatively 


because of the geometry of the sphere 


kro + = = Xr +o 


12 24a” 
= ka 6 = & qf ee e + .:- .) (50) 
24 24 


= kaé. 
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Then it immediately follows from Eqs. (9), (31), (35), (48), and (50) that 





ikaé —in/6 
,) € ix/ exp (— Trl €) 
V(P) = ” ~ 2e : i ix/3_-1 xs 2/3, -2° (51) 
0 n=0 € T T. (2a) v 


Once again we have a solution which we can compare with the classical solutions 
for electromagnetic radiation; and happily enough we shall find much that is similar. 
Indeed, van der Pol and Bremmer [21] in finding the Hertz potential excited by a vertical 
electric dipole at the surface of a spherical, homogeneous earth, have suggested an 
equation which is—even to the extent of defining 6 as in Eq. (3)—identical with Eq. (51). 
Furthermore, Eq. (49) is identical with their so-called tangent approximation to define 
the 7, . It is of some interest to note that Fock [22] in still another attack on this same 
problem agrees even more directly with our results. He derives a contour integral which 
is very similar in form to Eq. (34) and then deduces a series expansion which is again 


identical with Eq. (51). 
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—NOTES— 
A METHOD OF VARIATION FOR FLOW PROBLEMS—II* 


By A. R. Manwenn (University College, Swansea, England) 


Summary. The method of variation of reference [1] is developed afresh in a slightly 
different manner which enables the main principle used in [1] to be derived directly 
and also makes the actual calculations much simpler. It is shown how a variety of 
problems concerning aerofoils possessing minimal properties may be reduced to the 
solution of integro-differential equations which determine the mapping of the aerofoil 
onto a circular region. It is briefly indicated how the method may be extended to three 
dimensional flows 

1. Introduction. In [1] the author has given an elementary method of variation suit- 
type of extremal problems suggested by two-dimensional aerofoil 


able for treating 


theory. Briefly, the method is to make small elliptic bulges in the boundary and after 
calculating the changes of the flow functions to equate to zero the variation of the 
functional which is to be minimized, for the case of infinitely flat ellipses. This process 
Was justified by appealing to the fact that for such flat bulges the velocity changes 
as well as the geometrical changes were in effect infinitesimal. Such a physical argument 
being not quite convincing the author also tested the principle by showing that all 
small bulges in the hodograph plane gave equivalent results. It will however be seen 
that it is quite easy to prove that for infinitesimal velocity changes, although not for 
general perturbations of the boundary, the first variations may be equated to zero. 

Phe results of the method appear in a conveniently compact form viz. as integro- 


differential equations (sometimes just differential equations) from which the mapping 
of the aerofoil on the standard unit circle may be determined. In some problems it is 
necessary to make auxiliary restrictions on the aerofoil, such as limitations on the chord 
or the velocity in the field. 

In three dimensional fields the device of conformal mapping is of course not available 
but, analogous to the above equations for aerofoil problems, the method yields certain 
relations between geometrical and field properties on the boundary of the field. A dis- 
cussion of the determination of the field from such conditions is left for a future paper. 

2. The method of variation. In [1] the following principle was treated as physically 
obvious: 

If a functional of the geometry of a closed curve and the velocities of an associated 
hydrodynamie field is maximized by a certain curve, this functional is stationary for 
all variations in which both physical coordinates and the velocities are changed in- 
finitesimally. 

On the other hand small bulges giving rise to finite changes of velocity will change 
the functional by a small quantity of the first order. This situation may be illustrated 
by the following simple case. It is easily shown that the ratio of the area of an ellipse 
to the strength of the doublet, giving an equivalent disturbance in the stream at large 
distances, depends on the shape of the ellipse. Let it be admitted that there is some 
aerofoil problem in which the ratio of area to disturbance at infinity is stationary (c.f. §5). 


*Received Feb. 12, 1951. 
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Then by taking two small ellipses of the same area it would be possible to find an equa- 
tion corresponding to zero variation of area but a non-zero disturbance at infinity. 
Equally one could determine a variation corresponding to no disturbance but non-zero 
change of area. The ratio cannot be stationary for such variations. Further illustration 
may be found in the detailed calculations of [1] for small elliptic bulges. It is also in- 
teresting to note that, since velocity changes are finite, small elliptic bulges give some- 
thing more than a second variation and so the sign of the changes is rather strong evi- 
dence, although not proof, of a true maximum or minimum. 

The principle at the beginning of this paragraph will now be embodied in a simple 
lemma capable of direct proof. 

Main lemma. If y satisfies a second order partial differential equation of elliptic 
type together with the boundary conditions on a simple closed curve C (and suitable 
auxiliary conditions) whilst J is a functional of C, in the geometrical sense, and of the 
velocities at points in the field and on C, then, for variations giving infinitesimal velocity 
changes everywhere, 6J = 0 if J is a maximum or minimum. 

Proof. If a small perturbing stream function is added to y it is readily seen that 6y 
satifies a homogeneous linear partial differential equation and so, if é6y is a solution so 
is —é6y. Now the new boundary C’: y + é6y = 0 is to be derived from C by drawing a 


normal to C of length 
ss C 9 6 
in = -| ov / | = -( v) (2.1) 
on _}. Us/e 


Except at the isolated stagnation points the velocity qg, is not zero so the perturbation 
is infinitesimal provided 6y is made zero at such points. 

Geometrical changes may be expressed linearly in terms of 6n. For example, if A is 
the area enclosed by C whilst Z is the length of C and ds the line element 


64 = | énds, (2.2) 


6L = | K énds, 2.3) 


where K is the curvature of C. 

It is therefore clear that both the velocities at a given point and geometrical changes 
are linear in 6y. This is also true of velocities on the variable boundary C and for the 
present purpose this is more important. 

Thus, the variation of velocity on y = O may be found as the sum of : 

(1) 0/dn(6y) calculated at C, 

(ii) the change in the undisturbed field due to displacing C a distance 6n (i) is clearly 
linear in 6y and (ii) is linear in 6n and so again in dy. If (d°y/dn")¢ were zero then the 
contribution would be zero to the first order. In fact, for elliptic equations, this case is 


easily shown never to arise. 


If therefore J is a maximum or minimum the existence of a — 6/ for every +6] 
shows immediately that for such variations 6J = 0 as in the better known geometrical 


problems of the calculus of variations. Simple as the proof appears the result is not 
trivial since it implies such results as those of §5 [1] which are difficult to derive by 


direct calculation. 
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3. Laplace’s equation in two variables. Let z be the aerofoil plane, ¢ that of the unit 
circle and ¢ that of the strip. For non-circulating flow, U being the stream velocity, the 


complex potential is 


w= Ut = UE+ 7"). (3.1) 
A su:table perturbing function is 
bw = w?U Do a,t"’, (3.2) 
where y» is small and the a, are real so that J(éw) = 0 on the real axis. 


If the image in the ¢-plane of the disturbed boundary streamline C’ is given by the 
vector ¢ = e’ [I + 6R] it follows that 


op — #1(8) 
= Sin (3.3) 


where 

n(0) = >. a, sin né. 
Following the method of the previous section, the change of velocity is the sum 6,v + 62, 
where 


jv=p > na, sin n6/2 sin 86, (3.4) 


6,v/v = —dR(1 + a) (3.5) 


F dz dz 
a= Rs 43/2), 
dé d¢/. 
(3.5) being most conveniently found by writing 


dw _ | dw / dz 
dei |? le dé \ ee 


ec 


with 


a 


and expanding numerator and denominator. 
If C’ is mapped on the unit cirele | ¢’ | = 1 in a new plane the complex potential is 
known in two forms and by comparison of both it and the velocities on C’ found from 


the two methods it is quite easy to establish the relation. 


50(0’) 


1¢ 
log fs; | = A@ cot 0’ 6’)? 


dy’ 
where 
2sin 0/A@ = uw’ >. a, cos né’ and c=e. 
These relations together with Poissons’ formula for the function log (d¢/dg’) regular in 
¢’| > 1 and with the real part given on | ¢’ | = 1 determine the mapping of the new 
region on the old. For applications however 3.3, 3.4, 3.5 are sufficient and involve only 
one linear transformation of 7(@) whereas the explicit mapping requires three such trans- 


formations. 











108 NOTES [Vol. IX, No. 4 


4. Aerofoils minimizing a surface integral. Let A be the area contained in C a sym- 
metrical profile with axis along the stream and suppose 


/ = | F v) ds. 
where v is the velocity on the boundary. Then J TA~** is non-dimensional and it 
has been shown in [1] that for F(v) = kv with k constant J is a minimum for the circle. 


In fact J is not changed by conformal mapping between the aerofoil and plane of the 

circle but depends only on the potential function. The problem therefore reduces to 

making A as large as possible under conditions for which the circle is well known to 

be the result. It is natural to assume the existence of solutions for more general F and 

the method of variation shows that these must satisfy a certain integro-differential 

relation. From (3.5) since the equipotentials are normal to C it readily follows that 
6(ds) = 6.v'/v and so 

ds, = dé de fl + 6R(1 + a)! (4.1) 

de 

is the length of the element of the perturbed arc. Then 6/ [ F’(v)év ds + | F(v) 6(ds) 


becomes, after using (3.4), (3.5) and (4.1) 


é] = [ F’(v){ 6,2 - 6R(1 + a)}2sin 6d6 + ll F(v) 6RO + a) Zan fj 1.2) 
Again 
6A = i 6R = dé = Qu [ n\ 0) = ; 2 (4.3) 
4 dé J i 
so that if 6J = 0, which implies 
6] = KéA with K : 7 
the condition for a minimal solution may be written 
‘¢ H’(@)n(0) dé = [ F’(v) bv, dé, (4.4) 
where 
H = [ {fl + al[F’(v) — F(v)/v] + 2k sin 6 v’! dé. (4.5) 
Then 6,v, , the change in the plane of the circle is with v(@) = 7’(@) given by 
5,v, = Z [ oe an = ra dt (4.6) 


in which to make the operation on v possible it is supposed for example that v(/) possesses 
a derivative. 
Since 7 vanishes at the stagnation points at the ends of the aerofoil 4.4 becomes after 


an integration by parts 
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ee top dia = ; 5 
Pro] ; —= v(t) ar| do+ | H(t)v(t) dt = 0. (4.7) 
Jo Jo cost — cos 6 Jo 
If now v(¢) vanishes except in a small interval (4.7) gives 
“7 BI @)) «i 
3, FO) sin 6 ag +. HH) = 0, (4.8) 
Jo cost — cos 0 


where H is given by 4.5). 
In the case F = kv H simplifies greatly and (4.8) gives 


‘ : ( ‘ 
2k sin t id 1 — cost 2 
3 = { log ———i? =, 
v(t) mw dt | 1+ cost 7 sin ft 
Hence v« sin ¢ and | dz/d¢ | = constant showing that the solution must be a circle. The 
following approximations are suggested in other cases. 
(a) If the body is nearly circular a@ is small and the relation Eq. (4.8) is satisfied 
approximately by 
F(v) = v+ale" — 1) ve" = 2asin 8, 


where 6 is small. This gives one family of solutions. 
b) For thin flat aerofoils (i + @) is small except at the ends and if F(0) = 0 the 


term F — vF’ = 0(@) at the ends so that the equation reduces to 
"' 2k sin t 1 f" sin OF’ {v(6)} 

| ; dt + | : 1“! dg = 0. 
v Jo cost — cos 0 


5. Problems with restrictions. In this section a brief account will be given of the way 
in which the method can be modified to take account of certain restrictions peculiar to 
flow problems. For example, let A be the area of a symmetrical profile and D the strength 
of the equivalent doublet which represents the flow at large distances. As in §4 


? 


1 . f* (0) | dz |° * 
5. = q —_ ; — 1 . (e .1) 
on hg i sin 6 dé - : 


and since 6D is just the coefficient of ¢~' in the perturbing potential 3.2 


9 2 
sp = + | sin 6n(@) dé. 
T 


“0 


If no restriction is made the ratio D/A is stationary only if 


dz |* = i - dw 

; = Asin’ 6 or if 

d¢ dz 
is constant over the whole profile. This is possible only for a strip and to avoid this 
trivial result, where a finite area would have to be extended indefinitely along the 
stream direction, it is necessary to limit the aerofoil in this direction. The simplest 
restriction seems to be that the aerofoil is not to lie outside of two lines drawn per- 
pendicular to the stream. The ends then must lie on these lines and the method of varia- 
tion is not applicable to such parts of the profile. Over the rest of the profile dw dz is 
constant as before and the complete solution is that the aerofoils belong to the Ria- 
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bouchinsky constant velocity series [2]. In such cases the area must be regarded as 
fixed and the theory of isoperimetrical problems used to establish the constant velocity 
relation. For given straight lines, i.e. for a given maximum chord, the area completely 
determines the solution. 

In the last case the restriction was geometricaf but it is also possible to limit velocity 
instead of the chord. In the typical case there are constant velocity portions over the 
middle of the aerofoil and the problem is to find perturbations of the whole boundary 
which leave the velocity unchanged on these ares, say the intervals (1/2) + \ and 
(3x/2) + A for a symmetrical aerofoil. A simple method is to take a small bulge on the 
free are and then a general perturbation of the constant velocity arc. In the case quoted 
this is achieved by way of the transformation. 


t=¢+ 4 = gin a(s + b) (5.4) 


and a suitable disturbed potential is 








, ; I] 2 l = a,8 " a 
eer rE . f -e f + 1/¢ — 2 cos x T > sin d | wer 
where the a, are to be determined by the integral equation for 
v(u) = > na, sin nu 
r(u) — — AY __ 14 o()} [| Kew, Ole) at 
1/1 — sin’ X\ cos” u Jo 
= #8 ; : np" sin nu = Glu, x), 
= 
where 
. 2 = sin nu sin nt 
and 


8 = sin \/{cosx + ~/cos x — sin’ A}. 
In general successive approximation would be required on account of the presence of 
a(@) in (5.6) but if the aerofoil is flattish or if \ is small (5.6) gives v(u) directly to a good 
approximation. 

The author has carried out calculations for the problem of §3 with F 
quiring in addition that the maximum velocity is slightly less than twice the free stream 
value. In this case small constant velocity arcs appear and the method of variation 
gives the velocity over the rest of the profile. These solutions found by taking v(u) = 
G(u, x) in (5.6) tend smoothly into the circle as \ vanishes. It may also be noted that 
the result of imposing restrictions may be to give mixed boundary conditions. For 
example in the preceding problem if a curve of given area z is required to lie between 
y = +C where C < 1 the circle is replaced by a solution flattened so that y = +C 
over the middle whilst the method of variation gives the velocity over the ends. This 
means that over one part of the aerofoil the magnitude, and over the other, the inclina- 


v but re- 


tion of the velocity is given. 
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6. Minimal problems in potential theory. In dealing with aerofoils it is natural to use 
the conformal mapping of the region onto a circular one but this is only a device. As 
the following shows the true significance of the method is that it gives a functional 
relation between quantities on the boundary of the field. Moreover the linear trans- 
formations which arise in the method may be more compactly explained in terms of 
potential theory using Green’s Theorem for the original and perturbing potential c.f. [3]. 


Let 
i 7 0 Y {fe 
I= | r(Q, “ Se) dSé, (6.1) 
Jy on 
where = is an equipotential: ¢ = F and for simplicity it will be supposed that ¢ = 
O(1/r) at great distances whilst Q expresses the dependance of F on the space-coordinates. 


s 


If ¢ is varied the new equipotential L’ is given by drawing from ¥ a normal of length 


in = -| be / 2] ; (6.2) 


Now 6(dS) = Kén dS where K is the total curvature of = at Q and the first variation is 


61 = | FK indS+ | F, indS + | Foun in dS + | F,.(8¢), aS (6.3) 

suffixes denoting partial derivatives. 

The last term involves (é¢), but this may be expressed in terms of S itself by intro- 

ducing a new potential @ which takes the value of F,, on 2. 
Then, by Green’s identity, 


; , . 
F , (é¢), dS = [ = (*) dg dS = — [ #,¢, 5n dS, (6.4) 
Js ed on d 

according to (6.2). 


Applying the fundamental lemma of the calculus of variations to (6.3) the condition 
6] = 0 leads to 


KF + F, + Fo.em — Pe, = 9. (6.5 
As a simple illustration let C be the capacity of an isolated conductor having a charge 
which gives it a potential E above that at infinity. If V is the volume J = C/V" is 
non dimensional. 
Then 
6V = | ands’, (6.6) 


y 


and for the variations just considered 


5M 
sC = rz 
where 
1 f a . 
5M = — < S 
om hor I an (ap) 


is the change of charge on the conductor. 
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Hence 


s( l [ 0 , 1¢ 
6c = — ‘ — (dy) di 
trek Jy © On 
l f do. 
eS — 0—¢ ds 


trk Jy On 


oC = ~ ’ 5 I, (ae) én dss 6.7 


on 


Fron 6.6 6.7) sinc nis arbitrary the condition 61 0 1s 

dy 

const. 
on 
. ; ' A 

sot! he surface n st pe such that electric! \ is distributed uniformly over its surrace,. 
Thus 1 the solution ere not known in advance this result would lead directly to the 
sphe e = the solut 


Conclusion. The thod of variation of [1] has now been presented in a form in which 


the author hopes it may be Oo! practical utility. It ‘s a direct method for aerofoils and 
the numerical work in solving the equations vhich arise is not greater than in finding 
pp e flows for given boundaries: In some Cases much less. On the other hand it 
Is 1 expected that fo} smooth changes Irom true minimal shapes no large changes 
will ap] This 1s inde d an essential feature of all minimal solutions. 

TI thor is indebted to Dr. W. H. d Fuchs for explaining to him something ol 
the beautit il] work M. Schiffer |4] on minimal problems in conformal mapping It 
would ¢ t] such problems as arise 1n analysis are far deeper than those needed 
In iv heo 1 latter the condition O01 simple boundaries does not enter except 
tri « and since it is the essential condition mM analysis the present method may not 
be of anv use. It mmediately evident that the proof ot the main lemma of the paper 
exclua it regions here 7ward ariations would require a second Riemann sheet 

The author also thanks Dr. P M. Davidson for suggesting a variety of minimal 
prol lems in potent | theory and in particular the capacity problem of $6 to which he 
originally gave this solution using instead of Green’s identity for the potentials a phy- 
sical interpretation Ol the problem togethe1 with the method of virtual work 
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WALL EFFECTS IN CAVITY FLOW—II* 


By G. BIRKHOFF. M.. PLIESSET anp N. SIMMONS (Harvard University; Naval Ordnance Test Station, 
Pasadena; Ministry of Supply 


1. Introduction. In Part I of the present study,** the problems of flow about a 
cavitating body symmetrically placed in a channel or in a free jet have been solved in 
the case where the cavity extends to infinity downstream. The infinitely long cavity 
occurs, in each configuration, at one particular cavitation number which is a function 


of blockage ratio in the first case and is zero in the second. At greater values of the 


* denotes origin 
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FIG, - CASE A. 


cavitation number, the cavity is of finite extent and a different analysis is necessary. 
The solutions of the corresponding problems with finite cavities are given in the present 
Part. 

The configurations examined are again two-dimensional, this permitting the employ- 
ment of conformal transformation technique. The body is taken in the form of a finite 
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lamina perpendicular to the stream, so that the physical features of the flow may not 
be obscured by mathematical difficulties. Explicitly, the cases treated are 

A. The cavitating lamina in an infinite stream; 

B. The same in a channel of finite width; 

C. The same in a free jet of finite width. 

2. Case A. The lamina in an infinite stream. This case, where the liquid has no outer 
boundaries, is taken first to provide a standard for the other two. Solutions of this case 
have been given previously by Riabouchinsky [1] and by Fisher in an unpublished 
British Admiralty report, but the present treatment is much simpler than either. 

Take the density of the liquid as unity, the velocity at infinity as unity, and the 
width of the strip forming the body as 2 units, so as to avoid unnecessary symbols. 
This strip is disposed between the points (—1 + 7) in the z-plane (Fig. 1). The free 
boundaries, there shown in broken line, start from the edges of the strip and re-form 
downstream on a similar, conventional, solid strip, extending between the points (1 + 2). 
This device, which is due to Riabouchinsky, avoids the closure jets and turbulence that 
would otherwise have to be taken into account. In this way, the mathematical ad- 
vantages of a symmetrical problem are obtained merely by modification of the down- 
stream conditions, to which the flow around the cavitating body is known to be in- 
sensitive. That this is so, has been clearly demonstrated by Gilbarg, Rock and Za- 
rantonello,t who, in an as yet unpublished analysis of the similar problem with down- 
stream closure by a re-entrant jet, find for low and moderate cavitation numbers, cavity 
boundaries and drag coefficients virtually indistinguishable from those that result below. 

The flow being symmetrical about the x axis, consideration is restricted to the upper 
half z-plane. The corresponding regions in the W and ¢-planes are shown in Fig. 1, 
together with the auxiliary plane of u. Symbolism is as in Part I. 

Proceeding by Kirchhoff’s method for discontinuous flows, the transformation rela- 


tions are found: 


+ itanh Ben x | (1) 


— ztanh 6 tan u 


¢ = —(1 + Q)' E 


dz ; 1 — itanh B tan u |'” 
= §,(8) cos df a | : (2) 





du 1 + z¢tanh £6 tan u 
where 
B = 4 log (1+ Q) (3) 
and S,(8) can be evaluated in terms of standard Jacobian elliptic functions of modulus 
k = sech B as 


a a’ 
i(B) k ad - E a kK ” 


a) 


The integration of (2) between appropriate limits then yields the cavity half-length and 


half-width a: 


+Partial results were detailed previously by D. Gilbarg and H. H. Rock. Nav. Ord. Lab. Memo. 8718. 
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E — k?K 
1 — S.(8) = k’? + E’ ae k Pa? 
(5) 
Toe __ ki — KF) 
a ] = S;(8) = Lo + <a i $e ad 


The intrinsic equation of the cavity boundary, referred to its point of departure C as 
origin is found to be 





tan @ 
s= S( —_— ——. ——-—. |, 
“a (tanh” 6 + tan” 9)'”” | (6) 


In Cartesian parametric form, this is equivalent to 


x = S,(3)| cosh’ ab k, oF a” a) — sinh’ a-P(k, oe - 0) 
a _ 6 n\, (7) 


~ (tanh? B + tan? 6)'” 








S,(B) sinh* F mee | 
y = §S,(8) sin { "ia . 
Y , (tanh? 8 + tan? )'72 


where /, F are the standard elliptic integrals of the second and first kinds respectively. 
The drag coefficient of the lamina, based on unit velocity, is found after some re- 


duction to be 





Cp = (1 + Q)S,(8), (8) 
where 
S,(8) = pee oe (9) 
Referred to the velocity on the cavity boundary, (1 + Q)'”’, the drag coefficient is 
C, = S,(6@). (10) 
The functions S, , S. , S; , S, are all easily calculable and there is no difficulty in 


applying the foregoing solution in any numerically given case. The sub-class of cases 
in which Q is small is of especial interest: the functions then degenerate, giving the 
following simple results: 


Lo =e (A +2) +00), 


r+4\Q?"Q 
8 
(11) 
a 
Cp a gi + ® + 0%, 
C,= + 0(Q’). 


er 
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3. Case B. The lamina in a channel of finite width. 


.iG.a-tCase & 


Take the same arrangement as in 


Case A, but with the liquid confined between two parallel rigid walls, distant 2/ apart, 


with respect to which the body 


is symmetrically placed (Fig. 2). Again restricting con- 
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sideration to the upper half z-plane, the region in the W-plane is now an infinite strip: 
this, together with the ¢ and u-planes, is shown in Fig. 2. 
Proceeding as before, but with greater complexity due to the additional singularities, 


one finds the following transformation equations: 


,enu + ik’ sn u 


c= -—(1+ Q)” (13) 
dn u 
dz PAhksn A eno u — tk’ snu enu (14) 
du rl +Q)'"? 1—k’ sn’ Asn’ u 
2+ ¢ 
dn A TQ yy (15) 
Q) 
The constant / is not known ab initio, but must be determined at the end to conform 
} } OVC! 
Che integration of (14), between appropriate limits, then vields after reduction the 
0 pressions Tor the ecometrical characteristics: 
rT ~~ (J in oi eletteiei - a : aad - . 
. K'E(A) + (Eb K’)A| 4 cos. (ed A) — kK’ sn A, (16) 
2/ en A / en A 
h : 2h - 
[k°sn A —Z(A) de Al , (17) 
r(l + Q) ic 
Qh k’ ' " ? . ' 
73 fam A — tan (k’ se A)]. (18) 
rl+Q) kena 
he cavity shape is obtained in Cartesian parametric form as 
/ 2 T () - ] 0,(A +) 
= \/, sn Acd A — Z(A)iv + log " 
v I () 2 ~o (A —_ ss 
(19 
Q , . : vay , 
/ * itan (k’ se And v) — tan’ (k’ se A)], 
r1+Q 
here the parameter v runs from 0 to 2A. 
\oain. after considerable reduction, one finds for the drag coefficients 
h¢ : 
Cy = af +Q- Y tan '(k’ se A) |, 
(20) 


’ h Q 1 , | 
Cc, =2 — ‘ kh’ se 
I [ aa a ( c A) 


With these relations, the solution is formally complete in its barest essentials: if the 
complete flow pattern is desired, the 2, wu relation can be found without difficulty by 
integration of (14) and the velocity vector at any point is then given by (13). 
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The numerical solution for any given case involves some complication. Given Q and 
h, the value of k must be found from (16); successive approximation is the indicated 
method. At the same time am A is found from (15); the remaining results can then be 
evaluated. 
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PIG - CASE C. 


It is soon found, on trial, that solutions do not exist for all combinations of Q and h. 
For each value of Q, there is a limiting value of the blockage ratio 1/h that cannot be 
exceeded: this limiting value is given by 
(1) l+4Q0-(1+Q)'" 1 Q tan? Q 21) 

. = Sec Geer ease O° Se. (2 
 —_ 1+aQ r1+Q 2(1 + Q)'” 


It is easily verified that, in the limiting condition, the length of the cavity is infinite, so 


aa 


that the solution degenerates to that of Part I. Moreover, the liquid at infinity down- 
stream is on the point of cavitating. Hence the limitation is an inherent physical one. 
It bears some analogy with the choking phenomenon in a transonie wind-tunnel. 
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The limitation, at low cavitation numbers, is extremely harsh, e.g. at Q = 0.05, the 
blockage ratio cannot exceed about 1/1500 (cf. Part I, Sect. 2). Alternatively, for a 
blockage ratio of 0.05, the minimum cavitation number obtainable is 0.6. 

When numerical values are considered, it is found that for admissible solutions, the 
drag coefficients for any given cavitation number are virtually the same as in Case A: 
this is due to the very low blockage. The cavity tends to be larger than with infinite 
fluid, i.e., in effect, the cavitation number is decreased by the fixed boundaries, especially 
when conditions are nearly critical, but comparison of calculated cavity contours shows 
that this effect becomes appreciable only at points substantially downstream. 

4. Case C. The lamina in a free jet of finite width. Take the same arrangement as in 
Case A, but with the body symmetrically placed in a free jet whose width at infinity is 
2h units. Still restricting consideration to the upper half z-plane, the region in the W- 
plane is again an infinite strip and the transformation planes of ¢ and wu are as shown in 
Fig. 3. In these planes, it is necessary to take into account the points J, M at which 
the boundary stream-lines inflect. 

Taking account of these singularities and proceeding along the same lines as in the 


two previous cases, one finds the transformation relations 


— 1/2 Hite — |" 





dz _ __2hk___—*[{ Hu +i |” 
du x(1+Q)” | He —g) (23) 
where 
K 
B = — log (1 + Q) (24) 


and H, is the Jacobian theta-function constructed, like cd u, with modulus k. k in turn 
must be found from the complicated integral equation 
r(1 + Q)'” Fl A(iu+ ip)? _. si 
ete = | 8 sn 7u du. (25) 
2hk Jo LH(iu — 78) 
In this expression, the complex radical takes its first quadrant value. 
In terms of k and 8, the cavity dimensions are now found to be 


hk ** H(u + iB) + Hu — iB) 
-= —eaEe 7. ———s - i725 ’ 26 
L=704+07 J, (Ha + iB)Hu — iB} “ om (26) 





hk °K Hu + 78) — H(u — i8) ” 
~ = rap 72 ye OS eo ee : 2 
a oI, {Hu + iB) Hu — 7a}? 2% (27) 


The intrinsic equation of the cavity boundary is 


a ee 
~ (1 + Q)'” anv +kenv’ 

(28) 
tan 9 = 2 + 8) — HO — 8) 


i{H(v + 78) + He — 7p)}’ 


where the paremeter v runs from 0 to 2A. 
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The drag and lift coefficients reduce to 


’ Qhi , f’ HB + iu) — H(ib — iu) sn tu 
C; = 1+ Q)'” | Ue = = 3 = du, 
T Jy {Hsp + wH usb — tu)} 
29) 
Qhi * Hib + iu) — Hip — tu) sn iu 
CC; = | ——5 —— du. 
ri+Q)’ J, {Hap + wHasg — ru)} 2 


These relations comprise the solution of Case C. It is readily shown that, as / — 0, 
the solution degenerates to that of Case A. This however corresponds to very great 
values of h and is not of practical interest. In the general case, (24) and (25) must be 
solved simultaneously for # and 8 by a method of successive approximation. The re- 
maining results ean then, with some trouble, be evaluated. 

The case Q small, which is of the greatest practical interest, can be approximately 
solved in explicit terms. For one finds that this case corresponds to k — 1, so that the 
elliptic and theta-functions approach degenerate forms. Thus A is logarithmically large 


and 6/K small in comparison with unity. One develops the solution in powers of Q and 


retains terms of order Q. Then (24) and (25) become 
8 Q 
K r’ 
30) 
1+ Q) Q 
« \ B) t ~ \ 3 
h Tv 
where 
I l -+ sin ) 
S-(B ] Cos B sin B log . 
sin 6 
’ | 1 1 9 
S.(8 rsin 8 log see 8B + 2| f(tan 38) — J tan 50) | 
cos 8| f(sin 8) — /(—sin B 
in which 
Ha » 2) 
a (» | 
a tabulated function [2]. Hence, for given values of Q and /, 3 is readily determined 
The simplified forms of (26), (27). (28 29) sare. respect ely, 
2 l x : . : , o_ 
| asin 6 cos 8 log sin B + log tan =~, 5.3 
rl +Q@' LQ 2 
Zi - : : 
| cos B sin 8 log + sin B) |, 34) 
7 + @ LQ 
2h 
log cosh 
rl +Q 
0 sD 
) 
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Q. 
( 2h(1 + Q) [ cos B + E sin 6 log sec a|, 


T 
(36) 
2h OQ. 
( 1+Q) 2 [ —- cos B + * sin B log see a|. 
Phi foregoing general solution for () small is bound by the condition that 8 should 


not be small in comparison with Q. This merely implies an upper limit to the permissible 
vidth of jet and is no handieap in practice. Within the practical range of blockage ratios 


d cavitation numbers, the solution holds good. 


( 
t 


When Q is very small, the following first approximations may be used: 

: S;(8 

h ee 

2h . , 

sin 8, 
Q 
(37) 

2h ' 3 

a= ( — 060 8. 
Q 

Cp C, = Dil — cos 8). 
When Q — 0, those results become those for the infinite cavity discussed in Part [. 


It is not part of the present object to give detailed numerical results for application 
to arbitrary configurations: these it is hoped to present elsewhere. 
Acknowledgement is made to the Chief Scientist, British Ministry of Supply, for 


permission to publish Part TI of this paper. The views expressed in the paper are those 
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A NEW VARIATIONAL PRINCIPLE FOR ISENERGETIC FLOWS* 
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Rubinovy and the present author,’ it is shown that the variational 
itional flows of a compressible gas can be generalized to isenergetic 


tions to be varied are the stream function and the density distribution. 
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For such flows, L. Crocco has introduced a new stream function, which depends on the 
entropy. The advantage of this apparent complication is to make the velocity com- 
ponents directly expressible in terms of the partial derivatives of this new function. 
Thus a single differential equation containing only this stream function can be more 
explicitly obtained. In this paper, it will be shown that the same integral used in earlier 
variational principles yields Crocco’s equation when his stream function is being varied. 

Following Crocco, we shall refer all speeds to the maximum speed attainable in the 
field. The pressure-p and the density p are referred to suitable units consistent with this 
choice of typical speed. Then the isentropic acoustic speed c is 


(<) 


) 

c=y¥ P (1) 
p 

for an ideal gas with ratio of specific heats y. The condition of constant energy may be 


rewritten as 


2 a 2 
c= 7 LL = wy, (2) 


where w is the total speed. If the stagnation density for the stream-line of zero (reference) 
entropy is taken as reference, the density is given by 


oe® = (Y= yr. (3) 


where S is the entropy and R is the universal gas constant. 
Crocco’s stream function W is so defined that the velocity components u and v along 


the directions of increasing x and y are given by 
yul—w)’'” =¥,, 
(4) 
youl — wy)” = —-w¥,, 
where e = 0 is the two-dimensional case and e = 1 is the axially symmetrical case. In 


the latter case, as usual, z is taken along the axis of symmetry and y is in a perpendicular 
direction. The vorticity is 


e 2\ 1 4 1 ye ee l vy 
p27 - u=a¥ll— ww)” (W), Y) =- —— JP). 
w=? y= y w’) g(¥) yg, DR’ (WY) 


When wu and » are substituted from (4), this leads to Crocco’s equation for the stream 
| 


function V: 


( - wW. — 9, + (1 = Wn — et = yl — wy” 0 (8 -1). © 
Cc C Cc y Cc / 
We shall now show that 
bl = 6 [| (p + pw')y’ dx dy = 0 (6) 


with suitable boundary conditions, will lead to Eq. (5). By use of (1) and (2), the integral 


I becomes 
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: ¥ 
I = - || A(w”’) exp | -— [ g(¥) av |y dx dy (7) 
with 


A(w’) = {y¥ —-)N+04+ Dw} - w)””, 
where WV, corresponds to the streamline along which S = 0. In this integral, w is supposed 
to be expressed in terms of ¥, and W, through (4). In particular, w* can be expressed in 
terms of ¥2 + W? by 
y“ Bw) =Vv+V%3, Bw) =wi(l -— ww)” (8) 


The variation 6/ consists of two parts: (1) the direct variation of ¥, (2) the variation of 
WV through w*. The latter part can be easily transformed into variation of ¥ by (8). By 


noting that 
A'(w’)/B’(w’?) = yA — wy 


we obtain 


‘ 4534 ] i. 2 2y - ad 
= —— I A(w*) exp “Sa | g) dv |g) ov y° dx dy 


Y 1 Jy 


» nV 
+ [| (1 — w?) 4°) WL (dv), + ¥,(d¥),} exp | -— g¥) av |y dx dy. 
Jd i es Vo 


With boundary condition V,6¥ = 0, 6J = 0 leads to the equation 


: oi ei 2 .v 
“¢ ‘a — w)“O™ exp | -- Y_ |g) aw | #} 
Ox 7 — 4 y’ 
0 2 av Vv 
wy {a <r  O | -— gv) aw | “sh 
oY ie J o. y’ 


~ , + 
_ _A(w’) _ | | g) aw Joy. 
+= .> “Vo 


By direct differentiation and by use of (4) and (2), Eq. (5) is verified. 
For convenience of comparison, we record that 


dv S/R 


4. =e 


dy 


in the present dimensionless form, where y is the usual stream function. 
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ON MATRIX BOUNDARY VALUE PROBLEMS’ 


JULIAN ADEM ann MARCOS MOSHINSKY (Jnstitutos de Geofisica y Fisica, 
T 7] é sidad le VW érico 


By 


Introduction. In a recent publication’ a matrix type of boundary value problem was 
introduced in order to simplify the description of nuclear reactions. It appeared that 


this tvpe of boundary value problem could find applications in other branches of mathe- 


matical physics, and the purpose of the present note is to illustrate them. 

When we deal with vibrations of continuous media, with problems of heat flow ete 
we usually describe the state of the system in terms of a single function which depends 
on position as well as on the time. As an example, we may mention the lateral displace- 
ment of a vibrating string, or the temperature function in case of problems of heat flow. 

In many problems of 
below, the description of the state by a single function leads to boundary value problems 


vibration and heat conduction, of which examples will be given 


of great difficulty. It is possible though, in some cases, to divide the continuous medium 
into several regions, and with each region we can associate a function describing its 
state. These functions can be grouped together in the form of a column matrix or vector, 
which will then represent the state of the whole system. The mathematical problem we 


encounter then, is a matrix boundary value problem, which is, in general, much simpler 


than the one we would have to deal with in the usual formulation. 

In the present note, we shall discuss two examples of matrix boundary value prob- 
lems. The first one describing the flow of heat in a cross, illustrates the case where the 
interactions between the different regions appear through boundary conditions. The 


second one, dealing with the vibration of systems of plates with intermediate elastic 


media, illustrates the case where the interactions take place through the equations of 
shall obtain the eigenvalues and eigenmatrix functions corresponding to 


i 


motion. We 


this tvpe of problems, and with the help of them, give a formal solution for any initial 


conditions 

For the discussion of the self-adjoint properties of this type of boundary value 
problem, and the rigorous derivation of the series expansion theorems, we refer to other 
publications 

1. Flow of heat in a cross. We shall consider the problem of flow of heat in a cross 
Fig. 1a) whose four arms are of the same length /, and of square cross section of area a’, 


where a < The material of the cross will have a density p, conductivity « and specific 
heat «. The lateral sides of the cross will be coated in such a Wav that the outer con- 
ductivitv’ can be taken as zero, i.e. there is no radiation. 

It e tried to deal th this problem as a three-dimensional heat conduction problem 


in 2 region bounded by the surface of the cross, we would have a difficult boundary 
value problem which would not admit a simple solution. Taking into account though, 
ts us to assume that the temperature at all 


that the smallness of the cross section perm! 

points i t is the same e can describe the state of temperature in the cross in the 
following fashion vith each bar of the cross we associate its temperature Tunction 
6 her 2, 3, 4 indicates the bar in question, x represents the position of 
ne pon 1 re 0) I / as indicated in Fig. la., and / is the time. The 
temperature state of the whole cross is then described by the column matrix 
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6,(x, t) | 


| 


| 6.(x, t) 


-6,(x, t) 

The equation for the temperature in each bar will be the well known equation for 
heat flow in rods,’ and so the matrix representing the temperature in the cross satisfies 
the equation: 

pc(d0/dl) = x(0°0/dzx*). (2) 


The boundary conditions on the temperature at the free end points x = 0, can have 


Il 


abe 


| 


xp aaa 

















otk 
Nh 
‘ +h 
eH 

! 

| 

| 

! 
A 











C3 
_ 


iny of the usual forms;’ for simplicity we will assume that the end points will be main- 


tained at the constant temperature zero. We then have: 


6(0O, 4) = 0. >) 

We now consider the boundary conditions at the end points « = |. The intersection 

of the arms of the cross at x / gives rise to a cube of linear dimensions a illustrated 

in Fig. Ib. The smallness of the linear dimensions of the cube, allows us to assume that 

the temperature at all points of the cube can be taken as the same. The temperature 
at the end points of all the four arms will then be equal, and we have: 

A,(1, t) 6,(1, #) 6.(1, tf) = 6,(1, 0, 4) 


hich gives rise to three linearly independent boundary conditions. 
Finally, to obtain our last boundary condition we need to consider the total flow of 


eat into the cube. The flow of heat per unit time through the end sections of each bar 
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into the cube, is given by —xa°(d6,;/dx),-; . The net inflow of heat into the cube must 
be equal to the increment per unit time of the quantity of heat in the cube which is 
pca’(d0,/dt),.; . This quantity depends on a* and due to the smallness of the linear 
dimensions of the cube, it can be taken as of higher order than the net inflow of heat. 
The net flow of heat into the cube may then be assumed as ~ 0 and, as there is no 
radiation through the lateral sides, this leads to the boundary condition: 


4 
>> (00;/dx),-, = 0 (5) 
L=1 


The problem of flow of heat in a cross has now a complete mathematical description 
in terms of the equation (2) and the boundary conditions (3), (4), (5). Due to the sym- 




















Fic. 1b. 


metry of this particular problem, a simple linear transformation of the matrix @(z, ¢) 
can be found, which reduces (2)-(5) to four independent scalar problems. We shall 
discuss it though, as a matrix problem to illustrate the general procedure. 

If we now introduce, as usual, a solution of the form @(z, t) = @(x) exp (—At) where 
\ is an arbitrary real positive constant, we are led to a matrix boundary value problem 
in which 0(x) satisfies the ordinary linear equation: 

d°0/dx*) + (Ape/x)® = 0 (6) 

as well as the boundary conditions (3), (4), (5). 

To determine the eigenvalues and eigenmatrix functions of this problem, we first 


notice that from (6) and (3), 6(x) must have the form: 


0 . 1/2 ("7 
6(x) = Asin (Apc/x) “2, (7) 

where A is a constant column matrix of components A; , 7 = 1, 2, 3, 4. 
We now apply the boundary conditions (4), (5) to the solution (7) and we obtain 


the homogeneous system of linear equations: 


A, sin (Ape/x)'”7l — A, sin (Ape/x)'7l = 0; 1=2,3,4 
(8) 
4 


pe (Ape/x)'’? A; cos (Apc/ x)?" 1 = (). 


=] 
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The determinant of this system is: 4(Apc/«)'” sin® (Ape/x)'7l cos (Ape/x)'*7l, and 
the characteristic values for which this determinant vanishes are: 
22 
MWK 
= —; n=0,1,2,---. 9 
n 4pcl ’ , ? b ( ) 
We see from the determinant that when n is odd, the eigenvalue is non-degenerate, 
while when v is even, it is triply degenerate. The corresponding eigenmatrix functions 
are: 


for odd n: 





[2] 
| 
FE 
6,(z) = | (21)~*”? sin (nx2x/2l), (10a) 
1 
| 
Ly 
for even n: 
f 1] [ 0] 
| | | 
0| | 1 
0S (z) = 1”? sin (nwx/2l); 0. (2) = | \1-'? sin (nxx/2l) ~—(10b) 
“1 | 
| | 
| 
04 L—J4 
f 1 7 
| =i} 
0," (x) = | (21)~’”? sin (nra/2l). 
1 |} 
L—]J 


We define the scalar product of two matrix functions 6(x), (x) as: 


al I 4q 
(6,9 = | &@ue) ar = | b o.(a)y.(a) | dz (11) 

“0 “0 i=l 
where 6’(x) is the transposed form of the matrix $(x). It is then seen, that for even n, the 
eigenmatrix functions 0,°’(x), a = 1, 2, 3 were chosen so as to be mutually orthogonal, 


i.e. (05°, 0S”) = Oete. All the eigenmatrix functions are normalized (6, , 6,) = 1. Finally, 
the eigenmatrix functions corresponding to different eigenvalues, are orthogonal, as can 
be seen directly from (10), and also from general considerations of self adjointness given 
in another publication.* 

We assume now that the initial temperature distribution is given by a matrix function 
<(x) of class C"”, with sectionally continuous second derivative, which satisfies the 
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boundary conditions (3, 4, 5). The variation of temperature with time will then be 


represented *"* by a matrix function 
Kr. 2)-= Ss 6i4i8 x) exp [—(2n + 1)°axt/Apcl | 
(12) 
+ > 2 as? 05°) (x) exp [—n°a'xt/pcl’], 
where: 
a = (t, 6 Gs, = (#, 0, ). 


We see that in the present formulation, the problem of flow of heat in a cross admits 
a complete solution. 

2. Vibration of two circular plates with an intermediate elastic medium. Let us con- 
sider a system of two circular plates of radius R, clamped at the edges, and with an in- 
termediate elastic medium. We shall designate by p, the density, D, the flexural rigidity 
and a, , the thickness of the first plate, and p, , D, , a. will have the same meaning for 
the second plate. Finally, we denote by & the load per unit area of the plates necéssary 
to produce a unit compression in the elastic medium. 

The state of the vibrating system can then be described in terms of the normal 
displacements of the two plates u,(r, ¢, ¢) and u(r, ¢, t), in which 7, g stand for polar 
coordinates in the plane. As the force per unit area that the plates exert on each other, 
is proportional to the compression u, — U, of the elastic medium, we have that the 
equations of motion®’’ for the vibrating system are: 


p,a,(0°u,/d) + D,Ju + ku, — uw) = 0, 
(13) 
p2a(d°u,/dt) + D,Tu, + ku. — uw) = 0. 
We introduce, as usual, a solution of this system of differential equations, of the form: 
) COs m¢y| 
ul’, g, = ar exp (dwt) (14) 


= an 


in which, for simplicity in notation, we abstain from associating explicitly an index m 
with the two components column matrix u(r). We are then led to the matrix boundary 


value problem: 


dD, 0 | td 7 = u(r) wen — * k U(r) ‘ 
, ie = = § (15a) 
Yr ar ( r / ‘. 
0 D. Us(7) k W AsPp>o — k Uy(7) 
ul) = 0 (du/dr),-r = 0 (15b) 


This type of matrix boundary value problem differs from the previous one in so far 
as the interactions between the components take place through the equations of motion, 


and not through the boundary conditions. 


To find the characteristic frequencies and eigenmatrix functions of (15), we first 


consider the scalar boundary value problem: 
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(£? — r*)o(r) = 0; ~—-o(R) = O, (dv/dr),-2 = 0, (16) 


where £ is the operator [(1/r)(d/dr)r(d/dr) — (m*/r*)], and X is a parameter. The 


solution of (16) is well known, as it is the boundary value problem of a single circular 
plate.’ The function v(r) is given by AJ,,(\'?r) + BI,,(\'?r) and the eigenvalues 


\, , n = 1, 2, 3 are given by the transcendental equation: 
ee, ie ‘ . . a 
} Ja “r) (Mr) — I,,.(N?r) = J,(N”r) = @ (17a) 
L dr dr _— ; 
The roots 8, = Ai°R/x of this equation have been evaluated’ for several values of m. 
The corresponding eigenfunctions are: 
1B,1 J ,,.(wB,, 1rB,T i 
v.(r) = a 4.( ) _ AB) 1,(% a} (17b) 
R T,,, (B,,) R 
where A,, is an arbitrary constant. 
We propose now a solution for (15) of the form u,(r) = cv,(r) where c is a constant 


column matrix of components c, , ¢2 . The boundary conditions (15b) are immediately 


satisfied because of (16), (17). 
As £°v, = dAzv, , we see that the equation of motion (15a) is transformed into the 


algebraic linear equations 


([ Dias +k —k ap, 0 ( C; 
= ( (18) 
( —k D»x+k 0 Asp» Co 
To have non-trivial solutions of this system of equations, the determinant of the 
matrix must vanish, and this determines the two characteristic frequencies w,'', w, 
corresponding to each eigenvalue \, , which take the form 


( , ), 2 k ).2 tay 2 1/2) 1/2 
; (! £ Be 4. 5 Di:) - Es _ | \ , (19) 
“ a; py As Pe 1, Pp; A2Ppo 


1 


0! 
of 


Wn 
where 


a(n?) = (2a,p;)"(k + D,d2) — (2azp.)'(k + D.d;3) 


The corresponding eigenmatrices take from (18) the form 


k Dd, (A2p2 D.)y() 
ce’ = . c = 3 (20) 
—(a,p,/D,)y(\3) k/D. 
where 
v(A;) = —a(nz) + [a’(X;,) + (a, p,@2p2) ‘k*}' a 


The matrix boundary value problem (15) has now been solved, with the characteristic 
frequencies being w,'’, w,” and the corresponding eigenmatrix functions having the 
form: u!? (r) = ¢'”»,(r), uf? (r) = 2, (r). 

It is easily established, from the general equations (13) and the boundary conditions 
(15b), that two eigenmatrix functions u*, u corresponding to different characteristic 


( 
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frequencies, are orthogonal in the sense that: {% u*"Wu rdr = 0 where u*’ is the transposed 
of u* and w is the matrix 


a; py, 0 


0 As po 
, . : : . . . ~ Fi 
We can normalize the eigenmatrix function u in the sense that [*% u’Wu rdr = 1 by 
Ss J ° 
choosing the constant in (17b) appropriately. We would then have 


oR 
a 8) a)! 8 ' (6 
(u,,u; ) = u, Wu; rdr = dapdn , (21) 


where a, 6 = 1, 2, and n, 1 = 1, 2, 8, ; 

We have obtained an orthonormalized set of eigenmatrix functions corresponding to 
this vibration problem. With their help, we could represent the state of vibration for 
the two plates corresponding to any initial conditions. For example, let us assume that 
we had at ¢ = 0 a given displacement for our two plates, and that the initial velocity 
was zero. The initial displacement of the two plates could be represented as sum of 


terms of the form 
) cos me | 
ig ee 
sin mg 


For each term of this type, the solution of the vibration problem would be 


for m = 0, 1, 2, 


) cos me | 
wr, &) 
sin wal 


where u(r, ¢) has the form: 


u(r, 2) = >> {(e, wf? )ul?(r) cos wf? t + (2, ul? ul? (r) cos wi? t} (22) 


and 
sR 


(c,u, ) = | «'(r)Wu;"* (r)r dr. 


The generalization of the present developments to systems of more than two plates, 
as well as to other types of boundary conditions, and other forms for the plates, is 
straightforward. 

3. Conclusion. A general type of matrix boundary value problem, which includes 
both of the preceding cases, would be the one in which there are interactions between 
the components of the column matrix, both, in the differential equations and in the 
boundary conditions. If we introduce additional components into the column matrix, so 
as to reduce the system of differential equations to one of the first order, the matrix 
boundary value problems introduced in the present note, reduce to a form which has 
been extensively discussed by Birkhoff and Langer.” While, in their formulation, the 
self-adjoint nature of the present problems is obscured, their proofs concerning the exist- 
ence and properties of eigenvalues and eigenmatrix functions, as well as of expansion 
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theorems, apply to the problems discussed above. We are justified then in using formal 
expansion theorems, such as (12), (22) in the solution of problems of the above type. 

The authors are thankful to Ing. R. Monges Lépez, Director of the Instituto de 
Geofisica, for the encouragement he has given to the present research. This work was 
supported, in part, by the Instituto Nacional de la Investigacién Cientffica. 
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A NOTE ON A VECTOR FORMULA* 
By H. LOTTRUP KNUDSEN (The Royal Technical University of Denmark, Copenhagen) 


Of some vector formulas compiled in a recent paper’ the one discussed in the present 
note seems to be of general interest in field theory. 

1. Derivation of the vector formula. Let B(r) denote a vector function of the position 
vector r, satisfying sufficient continuity and differentiability conditions, and let A denote 
a closed surface and V the region of space bounded by this surface. Using conventional 
vector notation we may then state Gauss’ theorem in the following way 

| da-B= | wv-B. (1) 
“A “Vv 
Letting ¢(r) denote a scalar function and @(r) a dyade function, both possessing sufficient 
continuity and differentiability properties, we may derive the following equations from 


Gauss’ theorem 


| dag= | we, (2) 
JA Jy 
- a 2 
da-® = | dv 7 -®. (3) 
“A “Vv 


Substituting in equations (2) and (3) 


g = r-B, (4) 


® = 5B, (5) 


*Received May 9, 1951. 

1H. L. Knudsen, Nogle vektorformler og deres anvendelse, (Some vector formulas and their application), 
Fysisk Tidsskrift, Copenhagen, to be published. 

2M. Lagally, Vorlesungen wiber Vektor-Rechnung, dritte Auflage, Leipzig, 1944. 
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and using 


V -r= 3, (6) 
Vr=e, (7) 
where ¢ is the unit dyade, we obtain 
| dar-B= | d&V(-B) = | d&(Vr-B+ VB-r] 
vA “Vv “yy 
(8) 


= | dwe-B+ VB-r] = | dw (B+ VB-r}, 


| da:rB = | dv V7 -(rB) = dv [V-rB + r-VB] = | dv [3B + r- VB). (9) 
JV 


JA ’} “y 


Subtracting equation (8) from (9) we find 


} da-tB— | daB-r= | dv (2B +1r-VB— VB-r] (10) 


vA vA “V 


or introducing the unit dyade e 


[ da-[re — er]-B = | dv (2B + r-VB — VB-r]. (11) 
JA Jy 
We consider now the special case where B(r) is an irrotational field, i.e. we assume 
that 
VXB=0. — (12) 
When this condition is satisfied, the dyade VB is symmetrical. We have then 
r-V B= VB-r. (13) 


Introducing this equation in (11) we finally obtain the following equation for an irrota- 
tional field B(r) 


. 


; | da-|re — er]-B = | dv B. (14) 
2 Ja Jy 

By the theorem expressed through this equation the volume integral of an irrotational 
vector field over a region of space may be converted into a surface integral extended 
over the surface bounding this region. That this conversion can be carried out, is evident 
from the fact that an irrotational vector field B(r) may be expressed as the gradient of 
a certain scalar field, the potential g(r), 


B = V¢. (15) 


The conversion of the volume integral into a surface integral follows then directly 
from (2). However, the theorem expressed by equation (14) has the advantage that by 
using this equation we may express the surface integral without knowing the potential 
g(r) of B(r). 

The use of the theorem developed in this section will be illustrated by an example. 
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2. The force on a body in a central field of force. Let a body V, bounded by a closed 
surface A, be acted on by a central field of force. The center of the field is denoted by 
O and the force per unit volume of the body by f(r), where r is the position vector with 


O as origo. The force density f(r) may be expressed as 
f(r) = g(|r|)o, (16) 


where g(r) is a sealar function of the distance from the center of the field, and where 


o denotes a unit vector coparallel with r. 


The force F, with which the field of force acts on the body V, is expressed by 
F = / fav. (17) 
Jy 


A central field of foree being irrotational, we may convert the volume integral in this 
expression into a surface integral, extended over the boundary A of V by using the 
equation (14), derived in the last section. We find 


F =< i. de-tes > até. (18) 


Letting n denote the outward unit normal to A and denoting the scalar surface element 
by da, so that da n da, we may rewrite (18) as 
iF J ‘ 
F=>5 | [o cos (n, 0) — nj |r| g(\1r|) da. (19) 


. 


“vA 
The expression in the square bracket in this equation has the simple geometrical meaning 
demonstrated in Fig. 1. 


n @ cosi,9)-n 


Fia. 1. 


The application of the above developed expression (19) for the force on a body in 
a central field of force will be illustrated by two examples. 


The force inversely proportional to the distance. Let us first consider a central field of 
force, in which the force is directed towards the center and is inversely proportional to 
the distance from the center. For such a field of force the function g(| r |) is expressed by 


gir) = —Kir}7 (20) 


| 


where A is a constant. By substituting this function in (19) we obtain the following 
expression for the force on the body in question 


K f K f 
F= - : | [o cos (n, e) — n| da = “2 | © cos (n, 9) da (21) 








434 NOTES [Vol. IX, No. 4 


as we have 


| n da = 0, (22) 


“A 
this integral expressing the vector areal of the closed surface A. 

As an example of the application of the formula (21) we shall calculate the force 
on a sphere with center P and radius RF in a central field of force of the type discussed 
here, supposing that the center O of the field is situated on the surface of the sphere as 
shown in Fig. 2. From the symmetry it follows that the resulting force will be parallel 


n 








R da 


Fig. 2. 


with the radius PO to the center of the field; in computing the force we therefore only 
need to retain the component in this direction of each of the differential contributions 
to the surface integral. Using the symbols shown in Fig. 2 we find from (21) the following 


expression for the force F 


> pe 
F = 4 | t cos ; cos s 2k sin 0R dé = —KzrR’t, (23) 
where the first factor cos 6/2 stands for cos (@, t), the second factor cos 6/2 for cos (n, @), 
and 2rR sin 6 Rdé for da. In this equation t denotes a unit vector coparallel with OP. 
Through direct calaculation of the volume integral (17) we find, by dividing the sphere 
in conical shells with their apex at the center O of the field and introducing the angle 
a = 6/2 as an integration variable, 


ax/2 a2Rcosa 


| | tKe™' COS @ é da 2ré sin a dé = — Krk’t, (24) 


where cos a@ stands for cos (9, t), fda 27€ sin a dé for dv. 
In the example discussed here a double integral has to be calculated for finding the 
force as a volume integral, whereas by using the vector formula (14) we get the force 
expressed as a single integral. 

The Force Inversely Proportional to the Square of the Distance. For acentral field of force 
in which the force is inversely proportional to the square of the distance from the center 
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of the field and directed towards this center, the scalar function g(| r |) defined in (15) 
is expressed by 
girl) = —K |r|", (25) 


where K is a constant. The force F, with which this field acts on a body V, bounded by 
a closed surface A, may be found by substituting (25) in (19). We hereby find 
F = _# [ lo cos (n, 0) — n] |r|" da (26) 
7%" : ; 
For demonstrating the application of this formula we consider a sphere with center 
P and radius RF in a central field of force of the type discussed here. The center O of 
the field is assumed to be situated on the surface of the sphere. The previously used 
Fig. 2 may also be applied as an illustration of the present problem. For reasons of 
symmetry the force will be parallel with the radius PO to the center of the field; in 
computing the surface integral we therefore retain only the component in this direction 
of the differential contributions. We find from (26) 


> ef -1 
F = -3 | | cos 2 cos 5 — cos o|| an cos 4 2rR sin 6R dé = -2 KRt, (27) 
where the first factor cos 6/2 stands for cos (, t), the second factor cos 6/2 for cos 
(n, 0), cos 6 for cos (n, t), [2R cos 6/2] "' for |r |~' and 2xR sin 6 Rdé for da. 

From the theory for Newtonian potentials® it is known that the force F from the 
field of force discussed here may be computed by assuming that the mass of the sphere 
is concentrated to its center. As the volume of the sphere is 41/3 R*® and the distance 
between the center of the sphere and the center of the field R, we therefore find the 


force F expressed as 


F = —(KR (# Rt. (28) 


This expression is seen to be identical with the expression (27), found by using the 
surface integral method. The result from the potential theory used above was derived 
by carrying out a double integration. The calculation has consequently been simplified 
also in the present case by the use of the vector formula (14). 


O. D. Kellogg, Foundations of potential theory, Berlin, 1929. 


THE CIRCULAR PLATE WITH ECCENTRIC HOLE* 
By SAMUEL D. CONTE (Wayne University) 
1. Introduction. It is well known that the equation governing the small deflections 
of a thin uniform plate, assumed homogeneous and isotropic, is 
DAAw + P(z, y) = 0, (1) 


*Received March 20, 1951. 
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where P(x, y) is an arbitrary load per unit area acting normal to the upper surface, 
D = 2d°/3(1 — o°) is the flexural rigidity of the plate, d is the thickness of the plate, and A 
the Laplacian operator in two dimensions. All quantities appearing in this and in all 
succeeding equations are in dimensionless form. It is assumed that there are no forces 
acting on the edges of the plate other than the reactions due to the method of fixing of the 
plate. The state of stress and strain at every point of the plate can be completely de- 
termined when the deflection w of points originally on the middle plane of the plate is 
known. Complete solutions of the plate equation (1) under arbitrary normal loads and 
for various boundary conditions are known for thin rectangular [1], circular [2], or 
elliptic [3] plates. It is proposed in this paper to present a solution of this problem when 
the plate is bounded by two circles, one of which is interior to and eccentric to the other. 
2. Bipolar coordinates. The boundary conditions are considerably simplified when 
bipolar coordinates (a, 8) are introduced. The conformal transformation 
z=a2+ iy =c tanh } (a + 78) = c tanh 3, (2) 
maps the rectangular region [a, S a S a@, —7 S B S a] in the ¢-plane onto the region 
between two eccentric circles in the z-plane. As shown in Fig. 1, the boundary, of this 


region consists of the exterior circle a = a, and the interior circle a = a), with a > a, . 


ha 


ems 
eX 


B C.D» 








Fic. 1. 


The constant c is the pole of the system. The straight line segments AB and CD are given 
respectively by 8 = 0 and by 8 = z. Under the transformation (2) the Laplacian operator 


{3 a” 
= his —z fe 
A I (2, - 2), 3) 


becomes 


where (ch)? = (cosh a + cos 8)’ is the Jacobian of the transformation. If (hw) is taken 


as the dependent variable instead of w, the harmonic and the biharmonic equations 


become respectively, 


/ > > ; 

} oO Oo” , re) P 0 

Aw = scl € =) — 2sinha — 2sin B — + cosha — cos B/r(hw) = 0, 
Oa” + op” * Oa + 0p BY 


AAw 


Il 


a* 20° 3" a” a 
} 3 —. a ee Y Sige 1) cw = Q. 
da’ + 0a’ OB” + ap* 0am - Jom - mn 
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Thus the biharmonic equation becomes a linear partial differential equation with constant 
coefficients [4]. Separation of variables now yields as solutions the following series of 
biharmonic functions which have the period 27 in 6; 


hw = ie ¢,(a) cos nB + y,(a) sin nB, (4) 
where @¢ (a) = A, sinha + B,asinha + D, cosha + E,a cosh a, 
Yo(a) = 0, 
d;(a) = A, + Bia + D, cosh 2a + F; sinh 2a, 
Vila) = Ai + Bla + D cosh 2a + EF} sinh 2a, 
¢,(a) = A, cosh (n + la + B, cosh (n — l)a + D, sinh (n + La 
+ E, sinh (n — la, 
V,(a) = AJ cosh (n + lla + BL cosh (n — Ila + Di sinh (n + Da 
+ E/ sinh (n — Ila, 


and4A,,B,,D,,E,,A,, BL, D,, Bk. (n = 0,1, ---) are eight sets of arbitrary constants. 

3. Method of solution. It is now proposed to obtain a solution of the boundary value 
problem consisting of Eq. (1) and one of the usual boundary conditions. The most com- 
mon of these boundary conditions are those corresponding to clamped edges and to 
simply supported edges. For a clamped edge the boundary conditions in terms of (hw) are 


_ d(hw) 


— 0. (5) 


hw 
For a simply supported edge the boundary conditions are 
(hw) = 0, 


, 9? a” : 0 ‘ 0 
cosh of - +o (1 + oc) sinha = +(1+ o) sin B=) (6) 


Oa os” 


+ « cosha — cos ah i) = 0. 
A solution of the boundary value problem consisting of Eq. (1) and conditions (5) or (6) 
is sought in the form 
hw =wt+w, (7) 


where wy, is the series of biharmonic functions (4) and w,/h is a particular solution of the 
plate equation. If w, is expanded into a Fourier series, boundary conditions (5) or (6) 
give rise to a system of eight equations which can be solved for the eight sets of unknowns 
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which appear in w, . If, for example, a clamped plate is subjected to a uniform normal load 
p, a particular solution of the resulting boundary value problem is 


w A pc 
1-5 A=-— (8) 
and since the deflection must be an even function of 8, the biharmonic function is taken as 


UW, = Zz ¢,(a) cos np. (9) 


If w, is expanded into a Fourier series, the boundary conditions yield four sets of equa- 
tions for the coefficients A, , B, , D, , EF, which appear in ¢,(@). Explicit expressions for 
these coefficients are presented in the author’s thesis, and it is shown there that the series 
(9) with these coefficients converges absolutely and uniformly in @ and 8 over the domain 
between the eccentric circles. 

4. Particular solution for an arbitrary analytic load. Consider first a load function 
which is homogeneous of degree n in x and y. Expressed in bipolar coordinates such a load 


has the form 


Pla, B) = >> Qn smn —a SS . (10) 
k l 


A particular solution of the plate equation (1) for this load function is given by 


Ww : sinh” *“*? a sin‘*? 8 
= > A,,—— (11) 
h k=0 kh 7 


where the coefficients A,, are determined by the following system of equations 


2n—-—-k+2)\n—k4+ 1I(kK+ 2k + DAy 


+(in—-—-k+4)(n-—k+3)\(n —k+2)n —k+ DA, G-») (12) 
+ (hk + 4k + 3)(k + 20k + DAjuse) = Qu ; (kb = 0,1, --+ ,n). 
This system of (n + 1) equations can be solved uniquely for the (n + 1) unknowns 1,, 
(k = 0,1, --+ , n). These results for k = 2m, k = 2m + 1 are: 
\ _ Ll 1)’ (n — 2m + 2r + 2) oo+ (n — 2m + 3)(r + I) 
— (2m + 2) --+ (2m — 2r — 1) 


tr —2 


i _ (m Be 1)(r +. [n/2] ~~ mt 2h 
(r + 1)({n/2] + 2) n(2m 


r+metasa) (% + 2) --- (2 — 2m + 3) 
(2[n/2] + 4)(2m + 1)! 


(ry + 1)(2r + €) 
ae ae Qa(2tn/2)-2r)) 


"(n + 2) «++ (2[n/2] — 2r + 1) 
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ety (n — 2m + 2r + 1) -:- (n — 2m + 2(r +:21) 








Animes) = 2 a (2m + 3) «++ (Qn — 2r) 
J, (m+ D0 + [a+ 1/2] — m+ Die 
¥ (((m + 1)/2] + I(r + 1) n(2m—2r—1) 
(14) 
a (n + 1) --- (n — 2m + 2)(m + 1) 





my r+m+([n/2]—€a" 
+ 2 CD (Im + 1/2] + D@m +3)! 


_ (r + 1)(2r + &)! 
(n + 1) --+ (2[(n + 1)/2] — 2r) 





: Qn(2t(n+1)/21-2r-1)* 


In these equations {[n/2], for example, means the largest integer in n/2, and the 
symbol ¢, has the following meaning: «,, = 0 if (n + m) is odd; ex = Lif (n + m) is even. 
All factors in these equations must appear in descending order. When this is not the case, 
as for example, when r = 0, the entire factor is to be replaced by unity. 

Consider now an analytic load function P(z, y). If its Taylor’s series is rearranged 
in homogeneous powers of degree n in x and y together and expressed in bipolar coordi- 


nates, it has the form 





a. sinh"~* a sin* B — 
Pla, 8B) = >> > @u +" (15) 


n=0 k=0 


where @,, are known constants. By means of the method shown above a particular 


solution of the plate equation (1) is given formally by 


, eo a ee 
ud = > > ne sinh" * ae “8B. (16) 
The coefficients A,, for each n and k can be determined by means of Eqs. (13), (14). 

5. Special loads. When the load function has certain specialized forms, the peculiar 
properties of the operator A in bipolar coordinates may be profitably exploited to yield 
comparatively simple solutions of the boundary value problems in thin plate theory. 


It may be proved by induction or by direct calculation that if u, = sinh na/h", or if 
u cosh na/h”", then 
2n° 
a, = — is 
e 
and if v, = cos nB/h" or if v, = sin nB/h", then 
—2n’° 
ie = 6 « 
c 


Thus, if the load function is any finite or infinite sum of the functions wu, or v, , this 
stepping-down property of the operator A may be utilized to obtain a particular solution 
of the plate equation. A complete solution satisfying boundary conditions (5) or (6) may 
then be obtained by following the procedure given in Section 3. 

As another special case let the load function be 


P(a, 8B) = kh’ sinh pa cos 9B, (17) 
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where / is a known constant and p, q are integers but p # (q¢ + 1). If the plate is clamped 
at the edges, a particular solution of (1) is readily obtained in the form 


w,=A sinh pa cos gp, 


where 


D (p =e a 2¢° _ 2p" + | 
The boundary conditions may be satisfied by taking only those biharmonic terms involv- 
ing cos g8. Hence a complete closed form solution for a clamped plate under the load (17) 


is given by 
hw = A sinh pa cos g8 + ¢,(a@) cos gB, 


where @, = A, cosh (¢ + l)a + B, cosh (g — 1)a + D, sinh (¢ + la t+ E, sinh (¢ — la. 
The boundary conditions (5) yield four equations which determine uniquely the four 
constants A,, B,, D, , E, . Indeed, a solution in closed form can always be obtained if the 
normal load is any finite sum of terms of the following forms: h* cosh pa cos g8, h* cosh pa 
sin 98, h* sinh pa cos g8, h®* sinh pa sin gB. If p = (q + 1) in any of these loads the function 
is itself biharmonic. If, for example, 

P(a, 8) = kh’ sinh (¢ + Lae sin gf, (18) 
a particular solution of the plate equation which is periodic in 6 is given by 

a = Aa cosh (q rT l)a sin qe 

with 
k | 
D 24q(q + 1I)(¢q+ 4) 
A closed form solution involving only four of the biharmonic functions can again be 
obtained as in the previous case. We note, however, that if the plate is simply supported 


at the edges, no finite number of biharmonic terms will suffice to yield a solution of the 
resulting boundary value problem under loads of the form (17) or (18). 
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AN EXPANSION OF FREQUENCY DETERMINANTS WITH APPLICATION TO 

THE NORMAL FREQUENCIES OF A SPRING MOUNTED RIGID BODY 
(RESILIENT FOUNDATION )* 

By F. K. G. ODQVIST (Royal Institute of Technology, Stockholm) 


1. Introduction. The free small oscillations of a conservative system of n masses Q, 


with the coordinates g,; and velocities q; are governed by Lagrange’s equations of the type 
d & r" 


IW ; 
: — = 0, i= 1, ,n, (1.1) 
dt 09, 0g; 
where the kinetic energy is 
T=4)2 O24, 
and the potential energy is 


W=3 Do kaga - 


In the sequel we shall limit the treatment to the case where the masses Q; and the spring 
constants k;; are constants in space and time. 
The solution of (1.1) may be put in the form 


q: = PP aad 


Z hor? 6 (1.2) 
where A, is a constant amplitude, and w is the angular frequency the eigenvalues of which 
may be obtained from the equation 
—w +), Wi Win 
Ws —w +0», ie. 
= Q, (1.3) 
Wn Wr2 ; ——_ + Wan 
W here 
3 2 
we, 
Q; ‘ 


(1.4) 
There are also other mechanical systems leading to the same type of frequency equa- 


tion (1.3), e.g. the free vibrations of a rigid body, supported on elastic springs, if the 
coordinate system is properly chosen. 


In this case we have n < 6. 
2. Expansion of frequency determinants. A frequency or “‘secular’’ determinant of 

the form 

a, + 2, Qi2 Qin 

As; , a. + 2, Asn 
A(x) = 

Ant , Ane ’ 

*Received April 23, 1951. 
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may be developed in the series 





bad n 0, Aix | n 
A(z) = TT @,+2)+ 2] | TI” (@, +2) 
i i<k ap; ; 0 ‘ 
0, Ajt 5 Qj1 | 
n | n , 
+ Z, | ei P 0, Ay. | IT‘ aq. +2)+--- 
i<k<1 | | i 
Qi; ’ Aix ; 0 
0, Qi2 ; a 6 D2 Qin wi 
n Qo, ; 0, hia dan | 
+ 2] (a; + x)(a; + 2) 
s<i | 
ani 5 Ano 0 
0, Ayo a i 
ai Wer 5 0, as. 
+ del (a; + 2) 
Ant ’ On2 0 
| 
0, Qi2 ; Qin | 
Ao, , 0, Aan | 
+ ? (2.1) 
Ani 5 oe see 0 


which may be obtained by induction or by considering A as function of n variables 
y; = a; + x2 and using McLaurins development for this function of n variables.* 

In (2.1) upper indices on a determinant indicate that the minor, obtained by suppress- 
ing the corresponding rows and columns, is to be taken; upper indices on a product sign 
indicate that the corresponding factors in the product are to be replaced by unity. It is of 
interest to note that products of the type 


n 


II? @: + 9, 


a 
that is products of degree n — 1 in 2, do not occur in (2.1). 


*The latter method was suggested to the author by Mr. L. E. Zachrisson. 











1952] F. K. G. ODQVIST 443 


The expansion (2.1) may be used with advantage to approximate the roots of the 
determinantal equation (1.3). 

3. The roots of the equation A(z) = 0. Under the assumption that the diagonal 
terms a, , -:- a, of the matrix are all different and absolutely greater than the terms 
a;; , (¢ # j) outside the diagonal*, it is easy to obtain approximate expressions for all the 
roots of the equation A(x) = 0. 

To this effect we multiply all the terms a,;; with a common factor \, so that the 
expansion (2.1) takes the form of a polynomial of n:th degree in \, where the first degree 
term is missing. 

Putting the r:th root x, of the equation A(x) = 0 in the form of a power series in \, 
and determining its coefficients by introducing this expression into the equation, we 





obtain 
* | 0, ix —_ 
| IT"? (@, - @,) 
i<k | a; 0 | ‘ ‘ 
t, = —-a,+d foal n + r(---), (r= 1,--- n) 
(r) 
a; ~ @, 
II" (@ — @) (3.1) 
This power series in \ will converge for sufficiently small | \ |, which means a condition 


for the magnitude of the matrix elements | a,; |. 

The condition, imposed at the outset, that all the diagonal terms a, be different and 
dominate the a,; is not essential. 

In fact, it is easily seen that in the simple cases n = 2 and 3 it will make no difference, 
if one of the diagonal terms be included among the small terms, e.g. if one of the terms a; 
be multiplied with \ just as the terms a;; above. In these cases similar, though somewhat 
modified expansions for the roots x, will be possible. On the other hand it may be seen 
in the case n = 3, that if two of the diagonal terms, e.g. a, and a; , be multiplied with A, 
but not a, , then an expansion of the type (3.1) will be possible for x, but not for x, and 
x; , unless certain conditions be imposed upon the matrix. 

We shall leave the question at this point and give an illustration of the treatment of a 
special case of practical importance in the next section. 

4. Application to spring mounted foundations. 

Spring mounted or resilient foundations for stationary machinery, prime movers etc. 
are built in a great variety of types and forms. In many cases, though, the design will be 
something like the one shown in Fig. 1, and the assumptions made below will hold, at 
least approximately. We shall consider the foundation as rigid and the springs to have 
spring constants according to Hooke’s law. In this case there are in general six degrees of 
freedom of the system. Comparatively few references treat the problem with such 
generality [1]-[4], [6]. 

If the coordinate system be chosen so that its axes coincide with the principal axes of 
inertia at the centre of gravity C.G. of the foundation, the frequency equation of the 
motion will have the form (1.3). We shall assume the axis C.G.-3 to be vertical—which 
is often true or nearly true—so that the axes C.G.-1 and C.G.-2 are to be found in a 
horizontal plane through C.G. Further we take q, , 92 , g3 to be translations of the founda- 


*This assumption corresponds to what is generally termed a “‘weak coupling” of the system. 
tNumbers in brackets refer to the bibliography at the end of the paper. 
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dation in directions of the coordinate axes and q,/r; , qs/T2 , 9e/1s to be rotations with 
respect to the same axes, where 7, , 72 , 7; are the corresponding radii of gyration of the 


foundation. 
Furthermore, we assume the springs to have the following properties, often fulfilled in 


practice: 
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Fia. 1. 


(a) all springs have vertical axes, e.g. are statically submitted to pure compression; 

(b) each spring has a spring constant in vertical direction much larger than the one 
in a horizontal direction (for parallel displacements of the end faces of the springs). This 
will hold true f.i. for the ordinary type of rubber springs, used in this connection [5], 
[6], but it will generally not be the case for ordinary helical steel springs. Fig. 2 shows 
the relation between horizontal and vertical spring constant k, : k, for different values 
of the relative compression £ for the ordinary type of cylindrical rubber springs, vul- 
canized to steel plates, according to [6]. 

(ec) All springs have over-all dimensions small as compared with the foundation. If 
h,, , ho, , hs, are the coordinates of the point of application of the v:th spring and its spring 
constants in direction of the axes correspondingly xk,, , xk, , k;, where, according to (b), 
x may be taken as a small dimension-less quantity and k,, , k., and k;, are all of the same 
order of magnitude, then the springs may be arranged by the designer so as to fulfill the 


conditions 


Ksyho, _ 0, > ksh, _ 0, Dm ksh ha, = 0, (4.1) 
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the summations being extended over all the springs. From fig. 2 it is seen that « for 
rubber mountings, used in practice, will be of the order 0.1 — 0.15 (shaded area). 


kp 
ky 


O03 
0,2 


0,1 


-O2 





Fic. 2. 


From (a), (b) and (ce) it follows that the potential energy takes the simple form 


W = } > (xk,,51, + kk, 53» + ks, 53»), 


where 


9: + Qshz,/T2 — Qeho,/Ts , 


61, 
bo, = Qo + deh,,/T3 — Qahs,/T; ’ 


635 = ds + Qaho,/T - Qshi,/T2 . 


The assumptions (a), (b), (c) make possible further considerable simplification of the 


frequency equation (1.3). Due to (ec) we may write 


2 2 —2 
Wi; = wh; + Kay; , 
where the upper or lower sign should be taken so as to make @,; real. Then w* and @;, ; are 
of the same order of magnitude, if different from zero at all. 
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Accordingly, we put (Q = total mass of foundation) 


2 —2 K 2 -2 K 
wn = to =—~ > kh, Wis = Kois = ~— >, kyh , 
Q ’ Qr, ¥ 


| 
“ 
s 


2 —2 K 2 —2 K 
a> ee = i k,,hz Woo = KWo2 = 7 Zz, k, 
Or, S am os = 


2 —2 K 2 —2 K 
Gag EF “Ky Ee Or, : Koyhs, , Woe = KWog = Or, z keh, ’ 
1 v 3 + 
1 K 
2 . Po. —2 
a A be ks, ’ QS Ee ee py Koyhs, . 
Q , Qr, ’ 


] ; 

2 2 —2 2 2 

Wu = oh + kau = Or a (ks,ho, + kk>,h3,), 
1 v 





(4.2) 
2 —2 K 2 —2 K 
eS ee 2 Zz, k2,h, hz, ) @51 = KW51 = 7A bo ky hs, , 
Qnr, SF Qr. SF 
2 —2 K 
Ge = ~“KHy = - p> k, ,hohs, , 
Qrr, SF 
2 *%2 2 l 2 2 
W55 = Ws + KWs55 = C “2 > (ksh, + «k,,h3,), 
Yr a 
2 —2 K 2 —2 K 
HY = "ey = % D Kine» ] Wo = KWeo = = z. Koyhi, ’ 
Q) 3 - Qr, » 
2 2 K rs 2 —2 K 
Wes = —Kia = — 2 koh, ,h:; wes = — Kees = — A— D hilahs, , 
4 4 Ors : V2r'"l el ’Sy 9 5 € Qrors : Vie e2rl hs 
: ae ‘ ki ha, + Kophi 
Weg = KWeg = Qr2 3 (K,,N2, + Corltiy)- 
3 v 
. . . 2 . 
All the rest of the quantities w;; vanish 
2 2 2 2 2 2 2 2 2 2 2 2 
Wi2 = Wig = Wig = Wor = Weg = Wes = W31 = W32 = Wig = War = Woo = Wes = VO. (4.3) 
for reasons of symmetry, and 
2 ‘ 2 2 2 2 
34 = W35 = W43 = Was = W53 = Wa = O (4.4) 


due to the relations (4.1). 
If the expressions (4.2) to (4.4) for w;; are introduced into (1.3) the determinant will 


take the following form: 
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A(—w) = 
2 - —2 = 
—w +k; , 0, 0, 0, KWi5 » — KW16 
| 
2 —2 - 
0, —w + kw , 0, — Kes ; 0, KWo6 
0, 0, —w +w33 , 0, 0, 0 
—2 2 2 —2 —2 
0, — KW@40 , 0, a +wi; + Kw 44 ’ 0, —™ K@46 
2 2 2 —2 —8 
KWe1 5 0, 0, 0, ae +w*, + KW55 5 —™ KW56 
—2 2 ail —2 2 —2 
— KWe1 ; KWeo , 0, —a oe = Cites» —w + Kwes 


(4.5) 


Due to the comparatively great number of small terms w;; (¢ = 1, 2, 6) in the diagonal, 
it is not possible to obtain expansions in power series in «x of the type (3.1) for the roots 
of the equat ion. 

If, on the other hand, we submit the matrix to the additional conditions 


° if 


®i6 = We = @oo = @e2 — 0, (4.6) 


which are equivalent to the conditions 


be k. Rs, = 0, pm koh, = 0, (4.7) 


such developments will hold. It is to be noted that in many practical cases the conditions 
(4.7) will be a consequence of the two first conditions (4.1), due to the properties of the 
spring elements, such as disclosed by fig. 2. 

We may now put 


> 


2 2 2 —2 2 
oH = Kai, + Kin + ‘iota Wo = KWo. + K M2 + °°> 

= 2 i *2 8 2 
Os = W331], Ws = W44 + KW44 + K M4 + wit (4.8) 


2 2 —2 2 2 —2 2 
Ws = W35 + KWs5 + Kes tees, We = KWes + Ke + °°° 


Retaining the first term of the development (2.1), as applied to (4.5), untouched and 
ordering the rest of the terms according to ascending powers of x’, we may, after canceling 


the common factor w” — w33; , insert successively the expressions (4.8) in (2.1) and thus 
determine the unknown coefficients u, , --- , uw. Thus we finally get the simple solution 
—3 —2 —2 —2 
. = — Se b= WoaWa2 
Vibe 2 3 it *2 =? 
wis Was 
PF. ae xe 
46464 56/65 
ie = net cal ’ is ~~ e2 ’ (4.9) 
Was 55 


—2-—2 42 —2-—2 42 
_ 046064055 + WseWesWas 


x2, x2 
W44W55 
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5. Conclusions 

By using an expansion of the type (2.1) of the determinantal frequency equation (1.3), 
obtained for ordinary mechanical systems with n degrees of freedom, it is easy to obtain 
approximate expressions (3.1), supposing ‘‘weak coupling”’ and distinct eigenvalues. By 
examples it is shown, that similar developments of the eigenvalues as (3.1) will hold also 
in cases, where at least some of the terms in the principal diagonal of the matrix are of the 
same order as the coupling terms. 
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BOOK REVIEWS 


Electromagnetic problems of microwave theory. By H. Motz. Methuen & Co. Ltd., London, 
and John Wiley & Sons, Inc., New York, 1951. vii + 184 pp. $2.00. 


This book is a new Methuen Monograph dealing in Chapter one with a summary of general topics 
in microwave methods including velocity modulated tubes, travelling wave tubes, resonators, cavity 
magnetrons, methods of detection, wavemeters, standing wave meters, the Smith chart, and matching. 
Chapters two and three deal with a detailed discussion of velocity modulation and Klystron theory 
following the work of Webster. Chapter four is entitled Mode Selection in Cavity Magnetrons and 
includes the Fourier Analysis of the rotating waves and conditions for synchronism of the electrons with 
the rotating field components. Some discussion of mode stability is given together with an estimate of 
magnetron efficiency. 

Chapter five discusses the field relations in wave guides in orthogonal curvilinear systems. Applica- 
tion is made to guides of rectangular and circular cross section. The use of series expansions in terms of 
normal mode solutions is illustrated. The transmission line analogy is discussed briefly. 

Calculation of electromagnetic fields in cavities and guides of complex shape is discussed in some 
detail in chapter six where relaxation and finite difference methods are used. The methods are illustrated 
with resonator gaps. The determination of the field pattern for higher modes is also discussed briefly. The 
analytical treatment of corrugated wave guides, as done by Wilkinshaw, is outlined very well in the space 
of about six pages. 

Chapter seven is concerned with the impedance of an antenna in a wave guide. The problem of 
coupling into wave guides by means of straight wires or loops, as treated by the Toronto group, is pre- 
sented in a very satisfactory way. The theory of discontinuities in wave guides as worked out by Bethe 
and Schwinger is outlined and applied to transverse windows, changes in cross section, and bends. The 
methods are applicable only to thin windows and sudden changes in cross section but a numerical method, 
applicable to any type of discontinuity, is outlined. 

A large amount of useful material is packed into the 180 pages that make up this monograph. The 
writer did a very good job of discussing the problems which he includes. A person actively working in 
microwave theory will probably find this book useful as.well as the person looking for general information. 


Roun TRUELL 
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Ordinary non-linear differential equations in engineering and physical sciences. By N. W. 
McLachlan. Oxford at the Clarendon Press, 1950. vi + 201 pp. $4.25. 


The contents of this book are of essentially two categories. Following a general introduction, the 
second chapter deals with examples of non-linear equations which can be integrated explicitly. The third 
chapter is concerned with the particular situation where the explicit integration requires the use of elliptic 
functions. The remainder of the book deals with the second category, i.e., the non-linear equations which 
arise in vibration problems. 

In chapter four, the Van der Pol equation with small parameter, and the non-linear restoring force 
problem, are treated in detail. Emphasis is given to the physical phenomena involved as well as to the 
mathematical solutions. The relaxation oscillations of the Van der Pol equation are discussed very briefly. 
The method of slowly varying amplitude and phase and the equivalent linear equation are discussed in 
chapters five and six. Chapter seven is devoted to equations with periodic coefficients and chapter eight 
to numerical and graphical methods of solutions. Almost every equation in the book is illustrated by a 
physical problem so that the book should be of considerable interest to engineers. 


G. F. CARRIER 


Integraltafel. Zweiter Teil: Bestimmte Integrale. By Wolfgang Grébner and Nikolaus 
Hofreiter. Springer-Verlag, Wien and Innsbruck, 1950. vi + 551 pp. $5.80. 


This very complete table of definite integrals is a companion volume to Integraltafel (Unbestimmte 
Integrale) of indefinite integrals by the same authors, published in 1949 (reviewed QAM, January 1950). 
This book contains approximately two thousand definite integrals including some which are not found 
in Bierens de Haan. The symbols are large and clear and the general make-up of the table is very satis- 
factory. The only change that this reviewer would ask for is specification of the limits of integration in 
the table of contents. The table of contents is as follows: Symbols and Notation, Methods of Calculating 
Definite Integrals, General Integral Formulae. (1) Rational Integrands (detailed listing omitted). 
Orthogonal Polynomials, including Legendre, Hypergeometric, Tschebyscheff, Associated Legendre, 
Laguerre and Hermite Polynomials (detailed listing omitted). (2) Algebraic Irrational Integrands, 
including elliptic integrals of Legendre and Weierstrass Canonical forms (detailed listing omitted). 
(3) Elementary Transcendental Integrands. Integrands of the form R(exp dz, exp wr, °°), 
lexp (—sx)| f(x), R(x, exp Ax), R(x, exp f(x)), f(log x), log [g(x)], Euler Dilogarithm, f(x) (log x)" with 
f(x) rational, irrational and transcendental, f(x) log [g(x)], F[z, log f(x)], (exponential, logarithmic, sine 
and cosine integrals), f(sin x, cos x), f(sin az, cos bz, ---), F[z, sin az, cos bz], F[z, sin f(x), cos g(x), --*], 
Flexp az, sin bz, cos cx], F[x, exp az, sin bz, cos cx], F[z, exp f(x), sin g(x), cos h(x)], F[z, log f(z), 
sin g(x), cos h(x)], F[z, are sin x, arc cos z], Fiz, are tan z, are cotan x], R[exp Xz, sinh az, cosh bz], 
R[x, sinh ax, cosh bz], F[ f(x), sinh az, cosh bz]. Integrals of sinh 2, cosh z, tanh z, cotanh~ z. 
(4) Euler Integrals. Gamma function, Beta function, products of powers of linear expressions with 
general exponents, products of powers of binomial expressions with general exponents, products of 
powers of higher order expressions with general exponents. (5) Integrals of Cylinder functions. Bessel 
functions, Bessel functions with imaginary arguments, Integrals of the form F[z, Z,(zx)], F[x, exp z, log z, 
Z,(x)], F[x, sin 2, cos x, Z,(x)], F{z, Z,(x), Zu(x)). 


Roun TRUELL 




















